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ABSTRACT 

A  doubly  curved  element  for  a  shell  of  revolution  which  has  arbitrar- 
ily curved  meridians  is  developed  and  analyzed.    Meridional  curvature  is 
calculated  using  a  highly  accurate  polynomial  approximation.    The  dis- 
placement functions  selected  satisfy  interelement  compatibility  and 
contain    all  the  lower  modes  of  a  complete  set  of  straining  functions. 
Non-straining  modes  corresponding  to  rigid  body  motions  are  introduced 
into  the  final  stiffness  matrix  for  any  conceivable  rigid  body  motions. 

The  direct  stiffness  method  was  used  to  construct  a  stress  and 
strain  analysis  program  and  the  results  of  the  analysis  of  a  few  problems 
are  compared  to  classical  solutions  to  establish  the  integrity  of  the 
element. 
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I.     INTRODUCTION 

A.  GENERAL  DESCRIPTION 

The  concept  of  dividing  complex  structures  into  an  assemblage  of 
individual  structural  components  or  elements  is  axiomatic  to  the  struc- 
tural engineer.     However,  the  ability  to  execute  a  breaking  down  and 
reassembling  process  has  been  dependent  upon  the  availability  of    an- 
alytical tools.     Not  until  the  present  decade  were  developments  in  high- 
speed digital  computers  sufficient  to  make  possible  the  fundamental 
approach  to  problems  of  structural  analysis  known  as  the  finite  element 
method . 

The  finite  element  method  is  essentially  the  analysis  of  an  idealized 
structural  system  composed  of  a  finite  number  of  elements  interconnected 
at  a  finite  number  of  joints  or  nodal  points.     It  is  the  finite  characteristic 
of  structural  connectivity  which  distinguishes  the  finite  element  problem 
from  problems  of  continuum  mechanics  . 

B.  HISTORICAL  BACKGROUND 

The  finite  element  method  of  analysis  is  applicable  to  one,  two  and 
three  dimensional  structures.    Background  discussion  will  be  limited  to 
the  two  dimensional  case,  although  recently  there  has  been  considerable 
research  activity  on  three  dimensional  structures  . 

Various  shapes  of  finite  elements  have  been  employed  in  plane  stress 
and  plane  strain  analysis.    Early  research  was  focused  on  the  development 

of  flat  rectangular  and  triangular  finite  elements  .    An  abundance  of  papers, 
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technical  reports,  and  theses  were  produced;  a  complete  bibliography 
would  be  beyond  the  scope  of  this  paper.     It  is  sufficient  to  say  however, 
that  requirements  for  efficient  flat  plate  finite  elements  have  been  formu- 
lated and  tested.    Considerable  success  was  reported  in  Ref s .   1  and  2 
using  these  elements  to  approximate  shells  of  free  form. 

The  development  of  curved  finite  elements ,  a  considerably  more 
difficult  problem,  has  progressed  at  a  much  slower  pace.    Some  early 
success  was  achieved  in  the  development  of  truncated  conical  ring 
elements.    However,  problem  formulation  using  these  elements  has  not 
been  generalized.    Truncated  conical  elements  have  produced  excellent 
results  for  particular  problems  but  have  also  resulted  in  unsatisfactory 
analysis  of  similar  problems. 

A  ring  element  with  arbitrarily  curved  meridians  developed  by  R.  E. 
Jones  and  D.  R.  Strome  [Ref.  3]  produced  acceptable  results  for  the  case 
of  shells  of  revolution  under  pressure  loads. 

In  1968,  G.  Cantin  [Ref.  4]  successfully  developed  a  curved  finite 
element  for  cylindrical  shells.    This  curved  element  was  developed  using 
6  nodal  coordinates  per  joint  (uit  vs ,  wlf  w,t  ,    Wj-. ,  w^t-)  and  implicitly 
contained  a  complete  description  of  rigid  body  motions  .    This  element  was 
tested  and  found  to  produce  excellent  results  for  a  broad  class  of  cylin- 
drical shell  and  flat  plate  problems  for  both  membrane  and  bending 
behavior.     (The  flat  element  is  obtained  by  letting  R-*00.) 


C .     OBJECTIVE 

The  objective  of  the  author's  research  is  to  develop  the  stiffness 
matrix  of  a  four  nodal  point  finite  element  for  a  shell  of  revolution  which 
has  arbitrarily  curved  meridians ,  hereafter  to  be  called  a  curvilinear  shell 
finite  element. 


II.    THE  FINITE  ELEMENT  METHOD 

A.  GENERAL 

Finite  element  analysis  of  an  elastic  continuum  logically  separates 
into  four  phases:    (1)  structural  idealization,   (2)  evaluation  of  the  element 
stiffness  properties,   (3)  the  evaluation  of  the  force -displacement  be- 
havior of  the  overall  element  assemblage,  and  (4)     the     evaluation  of 
element  stress  and  strain. 

B.  STRUCTURAL  IDEALIZATION 

Structural  idealization  is  the  division  of  the  actual  elastic  continuum 
into  an  assemblage  of  geometrically  compatible  structural  elements  which 
are  interconnected  at  a  discrete  number  of  nodal  points,  through  which 
the  forces  are  transmitted. 

The  discretization  of  the  continuum  is  a  physical  process,  subject 
only  to  conditions  of  geometric  compatibility,  rather  than  an  approxima- 
tion of  a  mathematical  nature.    The  idealized  structure,  then,  is  con- 
structed of  elements  which  possess  the  material  properties  of  the  actual 
continuum . 

A  basic  reduction  is  necessary  to  make  the  idealized  structure 

mathematically  tractable.    This  reduction,  based  on  an  a  priori  selection 

of  possible  deformations  of  the  elastic  continuum,  constrains  the 

elements  to  deform  in  only  a  certain  number  of  predictable  shapes.    These 

constraints  called  displacement  functions  uniquely  define  displacements 

and  the  state  of  strain  within  the  element  in  terms  of  nodal  displacements. 

These  functions  are  the  keystone  of  a  successful  idealization. 

10 


C.  EVALUATION  OF  ELEMENT  STIFFNESS  PROPERTIES 

The  evaluation  of  the  stiffness  properties  of  the  element  is  the  crit- 
ical phase  of  the  finite  element  method.    As  the  evaluation  of  the  stiffness 
properties  of  a  curvilinear  shell  finite  element  is  the  object  of  the 
author's  research,  procedures  will  be  discussed  in  following  sections. 
For  the  present,  it  is  sufficient  to  say  that  the  element  stiffness  matrix 
relates  the  nodal  forces  and  nodal  displacements  .    This  relation  takes 
the  form 

[F]  =   [k]  {6}  (2.1) 

where 

{F}  =  nodal  force  vector 

{6}  =  nodal  point  displacement  vector 

[k]  =  element  stiffness  matrix 

D.  EVALUATION  OF  THE  FORCE -DISPLACEMENT  BEHAVIOR  OF  THE 
IDEALIZED  STRUCTURE 

The  force-displacement  behavior  of  the  overall  element  assemblage 
is  obtained  by  the  systematic  superposition  of  individual  element  stiff- 
ness matrices  and  load  vectors.    This  technique,  called  the  "direct 
stiffness  method,"  is  well  documented  in  the  literature. 

E.  EVALUATION  OF  ELEMENT  STRESS  AND  STRAIN 

When  nodal  displacements  are  known, element  strain  is  directly 
obtained  using  classical  shell  theory  strain-displacement  relationships . 
Element  stress  is  then  calculated  using  constitutive  laws. 
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III.     CONSTRUCTION  OF  THE  STIFFNESS  MATRIX 

A.  GENERAL 

The  stiffness  matrix  for  a  particular  finite  elment  is  the  relationship 
between  a  set  of  generalized  displacements  and  the  corresponding 
generalized  forces. 

Standard  matrix  notation,  square  brackets  denoting  rectangular 
matrices,  curly  brackets  representing  column  matrices,  superscript  (-1) 
denoting  inversion,  and  superscript  T  denoting  transposition,  will  be 
used . 

B.  GEOMETRIC  DERIVATION 

The  geometric  description  of  the  curvilinear  shell  finite  element 
must  contain  sufficient  parameters  to  guarantee  that  an  assemblage  of 
elements  will  form  a  smooth  surface  whose  slope  is  continuous  and  whose 
dimensions  at  element  boundaries  match  those  of  the  actual  continuum. 
Furthermore,  the  geometric  description  must  be  sufficiently  general  to 
describe  shells  of  revolution  which  have  meridional  inflection  or  reverse 
curvature  as  well  as  shells  of  degenerate  geometry  such  as  cylindrical 
shells  and  flat  plates.    An  extension  of  the  descriptive  procedure  used 
in  Ref.  3  was  found  to  satisfy  the  above  requirements. 

Geometrical  symbols  for  the  curvilinear  shell  finite  element  are 
defined  in  Figure  1.    Nondimensional  meridional  arc  length  is  denoted 
by 
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NOTE:  ALL  DIMENSIONS   BASED 
ON    SHELL  MIDSURFACE 

Fig.  I    DEFINITION  OF  GEOMETRICAL  SYMBOLS 
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*         l  (3.1) 

and  nondimensional  parallel  arc  length  is  denoted  by 

Tj=    *  .  (3.2) 

Through  the  use  of  a  second  order  polynomial  in  4  for  the  cosine  of  the 
angle  0  , 

cos  0  =  c0  +  c-yi  +  c2^  (3.3) 

and  the  values  of 

r'x',    r^2)   ,    0  M  ,  0^\jl,c  (see  Figure  1) 
expressions  for  completely  defining  element  geometry  were  directly 
formulated.     The  radius  of  the  element  is  given  by 

r(4)    =  r(l)    +   £j"    cos  0d^  =  r(0)    +  Mc04  +  ic142+|c248) 

(3.4) 
and  the  principal  radii  of  curvature  are 


Kl       d0/d4  ^  +  2^4  K      ' 


R2  =      rffl  „  (3.6) 

rf  sin  0 


where 


c0  =  cos  0^)  (3.7) 

gx  =,2  [3  f1 rJ~ )  "  2  cos  0(0  -  cos  0<2)]  (3.8) 

( 2)  _     (  i  )    v 

c2  =  3  [-2  (j —* j  +  cos  0(x)    +  cos  0(2)  ]  .  (3.9) 


14 


It  is  worthy  of  note  that  all  values  required  for  these  expressions 
are  found  from  engineering  drawings  or  models. 

For  particular  cases,  such  as  conical  shells  and  cylindrical  shells 
where  geometry  is  defined  by  analytical  expressions,  geometric 
approximation  is  not  required. 

C.     DISPLACEMENT  FUNCTIONS' 
The  displacement  functions 

(  U  \ 

U(4,T?) 

o*(4  ,n)j 

must  satisfy  three  critical  requirements. 

1.  They  must  insure  that  continuity  of  the  structure  can  be  auto- 
matically maintained  across  common  boundaries  of  adjacent  elements. 

2.  They  must  insure  that  all  the  lower  modes  of  a  complete  set  of 
straining  modes  will  be  present  in  the  stiffness  matrix. 

3.  They  must  contain  an  exact  description  of  all  possible  rigid 
body  motions  of  the  element. 

It  can  be  shown  [Ref .  5]  that  satisfaction  of  the  first  two  requirements 
will  guarantee  that  the  strain  energy  of  the  "assemblage  of  elements" 
approximation  will  represent,  for  specified  loads/a  lower  bound    to  the 
strain  energy  of  the  actual  continuum.     Satisfaction  of  the  third 
requirement  is  necessary  to  expose  hidden  constraints  opposing  rigid 
body  motion. 
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1 .      Application  to  the  Curvilinear  Shell  Finite  Element 

The  polynomial  expansion: 
11(4,17)  =    a:  ir)      +  a24  +  a3r)         +  a4  (3.11u) 

v(t,,n)    =  a54rj       +  a64  +  a7  77  +  a8  (3.11v) 

w(4,tj)    =  a9fr73    +a1043Tj2  +  an43r?  +  al243  (3.11w) 

+  a134ST)3  +  a144sTf  +  a1542T]  +  a16  43 

+  a17^T]3  +al84T72  +  aig4r?  +  a2c  4 

or  more  compactly 

[u]  =    [P]     {aj  (3.12) 

satisfies  the  requirements  of  interelement  compatibility  and  completeness 
of  strain  field.    Requirement  three,  that  an  exact  description  of  all 
possible  rigid  body  motions  of  the  element  be  represented ,  is  not 
satisfied,  however,  as  will  be  shown  in  the  following  section. 

D.     RIGID  BODY  MODES 

It  has  been  established  in  Ref.  6,  that  conditions  of  equilibrium 
are  not  satisfied  unless  the  displacement  functions  contain  an  exact 
description  of  rigid  body  motions  of  the  element.    Moreover,  it  has  been 
established  in  Ref.  7  that  inclusion  of  rigid  body  motion  is  essential  and 
cannot  be  represented  by  independent  displacement  components.    Rigid 
body  motion,  then,  must  be  introduced  without  compromise  to  interelement 
compatibility  requirements  . 
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A  general  method  for  including  rigid  body  motion  is  developed  in 
Ref.  6. 

1 .      Application  to  the  Curvilinear  Shell  Finite  Element 

Consider  a  point  P,  specified  by  position  vector  p,  on  the 
element  in  the  system  of  reference  shown  in  Figure  2.    If  the  element  is 
submitted  to  general  translation  then, 


6T   =  6ii   +  6sj    +  63k 


(3.13) 


u 


v 


w 
V    J 


-  sin  0 
0 

COS   0 


'V 


6S\        (3.14) 
v63> 


in  the  system  of  reference  (x,  y,  z) .    The  corresponding  scalar  compo- 
nents of  the  displacement  field  are: 

cos  0     cos  9         cos  0    sin  9 

-  sin  9  cos  9 

sin  0     cos  0  sin  0     sin  9 

or  more  compactly 

M  "    [RTRANS]  £6<3 
When  submitted  to  rigid  body  rotation  of  small  amplitude 

J3    =     fti  +  0a  I  +  £3k  (3.15) 

the  displacement  vector  field  of  the  element  in  reference  system  (x,  y,  z) 

is 


SR  =  j8  x  p 


The  corresponding  scalar  expression  is 
/    \ 


(3.16) 


u 


v 


r/ 


-  sin  0  (z  cos  0+  r  sin  0)       cos  9  (z  cos  0+  r  sin  0)     0 

-  z  cos  9  -  z  sin  0  r 

-  sin  9  (z  sin  0-  r  cos  0)       cos  9  (z  sin  0-  r  cos  0)      0 

(3.17 
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N 


N 


Fig   2  DEFINITION  OF  COORDINATE  SYSTEMS 
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or  more  compactly 

[u}  =    [RROT]  i8i] 
Combining  expressions  (3.14)  and  (3.17)  gives: 
\ 


u 


v 


w 


R, 


TRANS 


or 


R 


W  =   [R]  £u  "  } 


R 


ROT 


(3.18) 


(3.19) 


Matrix  [R]  is  given  explicitly    in  Table  1. 


E.      CLASSICAL  SHELL  THEORY 

In  the  introduction  it  was  emphasized  that  the  difference  between 
finite  element  problems  and  problems  in  continuum  mechanics  is  based 
on  a  physical  approximation  only.    Classical  theories  for  shells  from 
continuum  mechanics  are,  therefore,  applicable  in  the  finite  element 
approach. 

The  consitutive  law  for  an  isotropic,  homogenous,  thin  elastic  shell 
can  be  written: 

Kx        K2        0        0        0        0 

K2       Ki        0        0        0        0 

0  0        K3      0        0        0 

0  0         0        Di      D2     0 

0  0         0        D2     Di     0 

0  0         0        0        0        D, 


'N     ' 


N 


V 


N; 


iy\  = 


) 


M 


M 


T? 


.  M, 


f.,t\ 


r\ 


\ 


in 


(3.20) 


v 


IV 


J 
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where 

Ki   =  T^    ;  Ks  =  v  Ki  ''  K3  =   £  U^tj  Kx 

Di  =  i?^-tl5v'ps  =   ^Di    ''   D3  =  i  (l~v)  Di 

or  more  compactly 

{N}  *    [E]   {c}  (3.21) 

Young's  modulus  E  and  Poisson's  ratio  v,  according  to  definition 
are  constants.    The  thickness  t,  however/  is  allowed  to  vary  in  the 
meridional  direction.    In  this  thesis,  shells  of  uniform  thickness  and 
shells  of  linearly  varying  thickness  are  considered.1 

Numerous  theories  exist  for  curvilinear  shells.    The  venerable 
theories  have  been  subject  to  years  of  critical  inspection;  their  merits 
and  limitations  are  well  known.    Several  theories  investigated  would 
have  been  acceptable.    Kraus'  formulation  [Ref.  8]  was  selected  for  the 
following  reasons: 

1.  The  theory  is  consistent.  That  is,  the  theory  holds  for  the 
limiting  cases  of  cylindrical  shell  and  flat  plate  elements.    This  quality 
is  of  basic  importance  since  mesh  size  reduction  to  the  extent  where  the 
element  is  nearly  flat  is  necessary  in  many  finite  element  problems. 

2.  The  theory  is  presented  in  a  form  that  can  be  directly  applied 
to  this  study. 

3  .    The  theory  allows  strain  free  modes  for  rigid  body  motions  . 


1  Any  thickness  which  varies  as  a  function  of  £  could  be  considered 
by  changing  only  one  card  in  the  computer  program. 
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The  strain  displacement  relations  are; 


/        \ 


r\ 


\ 

x      I 

x4 


r\ 


x. 


4r? 

V     J 


i  34 


cos 

0 

r 

1 

3 

1 


ax    3  77 

.X  .A 

COS    0 

r  R, 


1 
Ri 

sin  (b 
r 


1 3r_ 

r4    3^ 


1         3Ri 
XRi2       3  4 


1 


42      34: 


cr    3r] 


1 


a  rRj,       37] 


_L  A_ 
4r     34 


3r 

rT4" 


1 

a2r2 

32 
hrf 

cos 
r  I 

<t>      3 
34 

2 

r2<*4 

3r     3 
34      3t? 

2 

32 

f     \ 

u 


V     \ 


r  a  i.343  77 


w 


(3.22) 

or  more  compactly 

{€}=    [D]  {u}  (3.23) 

All  classical  theories  investigated  by  the  author  contained  inherent 
singularities.    Kraus'  formulation  was  no  exception.    It  can  be  seen  that 
several  terms  in  equation  (3.22)  become  infinite  whenever  r  is  equal  to 
zero.    This  theoretical  inadequacy  has  been  avoided  by  specifying  either 
a  small  hole  or  a  small  rigid  insert  with  appropriate  boundary  conditions, 
at  the  pole  point.    The  problem  solved  in  Chapter  5,  Section  B,  2 
demonstrates  this  method. 
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F.  NODAL  COORDINATES 

For  computational  convenience ,  physical  nodal  coordinates  are 
transformed  to  a  generalized  coordinate  system.    This  transformation 
facilitates  the  expression  of  displacements,  elastic  characteristics,  and 
strain  energy  in  a  system  free  from  geometric  complications  .    The 

following  set  of  nodal  coordinates  was  selected: 

T 
{uj  =    <  u1#  vt,  w1#  Wj,s    ,    w1/Q  ,    wlfSQ     >  i  =  1,2,3,4 

(3.24) 

where  s  and  6  are  the  physical  coordinates  of  any  point  on  the  element. 

Then  using  (3.1),   (3.2),  and  (3.12)  the  transformation  matrix  [t]  is 

constructed  such  that 

{uj  =    [T]  {at}  .  (3.25) 

The  transformation  matrix  [T]  is  nonsingular  and  is  easily 

inverted . 

Matrix  [T]  is  given  explicitly  in  Table  2  . 

G.  STIFFNESS  MATRIX 

The  strain  energy  in  the  element  is: 

Us   =  i  JJ  {N}T    {e}  r  0!  £d4dr?  (3.26) 

Substituting  for  stress   [n]  and  strain     [e]  one  gets: 

Us   =  Hui)T  [T-1f(jj[P*lT[E][P*]ro  i    B^O^LT"1]  {ut}     (3.27) 

where 

[P*]  =    [D]  [P]  (3.28) 
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The  stiffness  matrix  is  then: 


[k]  =  [r1] 


[P*][E]  [P*]roJ  hibr)  [r1  ] 


or 


[k]=  [t-1]7  [k]  [r1] 


(3.29) 

where  [T]and  [k]   are  both  (24x24)  matrices. 

Stiffness  matrix  [K]  may  contain  some  rigid  body  modes  depending 
on  problem  type.  However,  by  subjecting  [K]  to  eigenvalue  analysis, 
individual  rigid  body  modes,  if  present,  can  be  identified. 

When  [k]   does  not  contain  all  six  components  of  rigid  body  motion, 
it  must  be  modified  accordingly.    The  method  developed  in  Ref .  6  is 


applied  as  shown  below: 


KI-  [ 


R 


--R- 

u. 

I ]  J 


(3.30) 


where   {u.     }    =     (  &x#  6    ,  6^  fix,  j8    ,  ^    >  ,       [i]  is  a  (24x24) 

identity  matrix,  and  [R]  is  the  (2  4x6)   matrix  developed  in  Section  D. 
Matrix  [k]  is  then  expanded  to  include  all  six  rigid  body  modes: 


0" 

where    [  [R]T  [k]  [R] 


R 


F 
0 

Solving  for 


[K] 


[K]       I    [K]  [R] 

Lr]t[k]  !  [r]t[k]  [R] 

is  a  (6x6)  matrix. 

:     [k]  lr] 


(3.31) 


[R]T[K]    i    [R]T[K]  [R] 


u, 
J 


R 


(3.32) 
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{u.R}     =     [[R]T[K][R]  J'  [[K][R]]    {uj 


(3.33) 


However,  if  the  stiffness  matrix  [k]  already  contains  some  of  the  rigid 

r       t 

body  modes,    then  matrix         [r]     [k]  [R]        will  be  singular.    These 

modes  have  been  previously  identified  and  are  merely  suppressed  in 
|  u.      [   and  in  the  corresponding  column  of  [R]  .     Matrix      [r]     [k]  [R] 

then  becomes  nonsingular.     Solving  (3.33)  and  substituting  into  the  first 
(3,32)  leads  to  the  following  expression: 

{F  }  =  (  [K]  -    [  [K][R]][[R]T  [K]  [R]  J"'  [  [R]T  [K]  ]  )     {  u,  } 

(3.34) 
The  stiffness  matrix  modified  to  include  rigid  body  modes  is  then: 

[K*]=   [k]  -    [[K][R]]  [[R]T[K][R]j_1  [[R]T[K]" 

Stiffness  matrix    [K*]  is  strain  free  for  any  rigid  body  motion. 


H  .     DEFINITION  OF  MATRICES 

In  the  following  sections  of  this  thesis,  the  stiffness  matrix    [K*] 

r       T 

derived  above  is  referred  to  as  CSFE.    The  matrix       [R]     [k][R] 

developed  in  (3.31)  is  hereafter  referred  to  as  K22.     For  the  general  case 
where  all  six  rigid  body  modes  are  absent,  K22  is  a  (6x6)  matrix.    In  the 
case  of  a  cylindrical  shell  segment  where  four  rigid  body  modes  are 
absent,  K22  is  (4x4),  and  in  the  case  of  a  conical  segment  where  five 
rigid  modes  are  absent,  K22  is  (5x5). 
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I.       LOAD  VECTORS 

Loads  and  moments  acting  on  the  element  surface  must  be  repre- 
sented by  vectors  which  are  consistent  with  displacements  and  rotations 
in  the  generalized  coordinate  system. 

A  consistent  load  vector  then,  is  defined  as  a  vector  made  up  of 
fictitious  quantities  which  would  give,  after  an  inner  product  with  the 
nodal  point  displacement,  the  same  work  as  the  real  load. 

1 .      Concentrated  Load  Vectors 

For  concentrated  loads  of    Ft  ,  F    ,   F^,  acting  at  a  point  (s,  6) 

of  the  element,  the  work  done  is 
W  =  F,u  +  F    v   +  F^w 

or  more  compactly 

W  =  {F,}T     {u} 
Substituting  for    [u]  from    (3.12)  and  (3.24)  gives, 

W,{P.f  [P][t*]{«,} 

The  consistent  concentrated  load  vector  is  then: 

{cload}=[p]T{f'}  <3-35' 

2  .      Pressure  Load  Vectors 

The  work  done  by  a  pressure  load  p  is: 

W  =         p  w   d  A  (3.36) 

JA 

then  if 
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{q}  =   <  0      0       p  > 
direct  substitution  of  known  quantities  into  (3.36)  gives 


W  = 


J{q}T    [u]  lctrdidr)=  (jLct  JJ{q}T  [  p]  [l"1  J  r  d4drj)  {  Ul 


The  consistent  pressure  load  vector  is  then: 


{pload}=  *«ll[pT  Mr.m,? 


(3.37) 
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Table  1.    Transformation  Matrix  [r]  Between  Nodal  Coordinates 
and  Polynomial  Coefficients 


Rl,l* 

0 

Rl,3 

0 

Rl  ,5 

0 

0 

1 

0 

R2,4 

0 

R2,6 

R3     1 

/ 

0 

R3,3 

0 

R3,5 

0 

*Si 

0 

R4,3 

0 

R4,5 

0 

0 

R5,  2 

0 

R5/4 

0 

0 

0 

^,2 

0 

^,4 

0 

0 

R7,l 

0 

R7    3 
/ 

0 

R7,5 

0 

0 

1 

0 

0 

0 

Ra,6 

R9,l 

0 

^,3 

0 

^9    5 
/ 

0 

Rio  / 1 

0 

f*L0  ,3 

0 

^0,5 

0 

0 

Rll     2 
/ 

0 

Rll,4 

0 

0 

0 

R12,2 

0 

R12,4 

0 

0 

R13,l 

R13,  2 

R13  73 

R13#4 

Rl3/5 

0 

R14,l 

R14,2 

0 

0 

0 

R14,6 

R15,l 

R15,2 

R15,3 

R15,4 

RlB/r5 

0 

R16,l 

R16,2 

R16  ,3 

R16,4 

R16  ,5 

0 

R17,l 

R17,2 

0 

R17,4 

^17,5 

0 

R18/l 

Ria,  2 

0 

Rl8     4 

Rj.8  f  5 

0 

R19,l 

R19,  2 

R19,3 

R19     4 
/ 

^9,5 

0 

R20/l 

R2o     2 
/ 

0 

R2o,4 

R2o,5 

R3o,6 

R21'l 

R21,2 

R21,3 

R21,4 

R21,5 

0 

R22,l 

R22     2 
/ 

R22,3 

R22  ,4 

R22f5 

0 

R23,l 

R23,  2 

0 

R23,4 

R23,5 

0 

R24,i 

R24,  2 

0 

R24    4 
/ 

R24;5 

0 

*    Commas  in  this  table  are  used  as  separators,  no  partial  differentia- 
tion is  implied. 
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where: 

Ri  i 

/ 

Rl,3 

Rl    5 

R2,4 

R2    6 
/ 

R3,3 
R3,5 


•(I) 


R 


COS    0  Vi 

-sin  0  i1  ) 

Rl,l   z~   Rl,: 

-z 

r(l) 

sin  0  (x) 
cos  0  (x) 


R3   j  z  -  R3  3^) 
Ri/i   ci 


4.1 


R3    1      I 


R 


4    3 


■(1) 


R 


4    5 


R3     1         "r     R^     -i       Z    T    Rt      ■,      C  _ 


ll.l       ^0 


R, 


5    2 


R 


s  4 


R3  .1 
R3,5 


R 


R 


6     2 


Ixi. 
(l) 


r 

R4    5 
^4     =     r(lJ 

R7  l   =   cos  0 


(2) 


R7  3  =   -  sin 
R7    5    = 


U) 


R73r(2> 


^8    6 
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Rg.i    - 

R. 

R, 


=   sin  0  (2) 


^9,3 


cos  <h 


(a) 


v9    5 


(3) 


Rio  ti    ~ 

RlO  f  3     =" 
R10  ^5       = 

RU     2    = 


"     ^9    3  F 

R7  x  (cj.   +  2c2) 


R9 
1 

li 

(Ci 

+ 

2c2) 

4 

R9/ 

2 

1 

+   R*,; 

Re, 

,i 

(cx  +  2c2)r 


(2) 


r(2) 


Rn  4  =  cos  0 


(3) 


R12,2 

= 

^0,1 
r(a) 

Rl2    4 
/ 

= 

RlOf5 

R13,l 

= 

(2) 

cos  0  v  '  cos  a 

Rl3,2 

= 

(2I    . 

cos  0  v  ;  sin  a. 

R13  ,3 

= 

-  sin  0(2) 

Rl3     4 

/ 

= 

-  sin  0^  sin  a  r^2' 

R13,5 

= 

sin  0  *s'  cos  a  A2) 

R14,l 

= 

-  sin  a 

R14,2 

= 

cos  a 

R14,6 

= 

r(s) 

R1S     1 
/ 

= 

(2) 

sin  0V  '  cos  a 

R15    a 

= 

(2)     . 

sin  0V  '  sin  a 

R15  ,  =   cos  0(2) 


^15,3 
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^15 


R15 
R 


16 


R 


16 


R 


16 


R 


16 


R 


16 


R 


17 


R 


17 


R 


17 


R 


17 


R 


18 


R 


19 


R 


19 


R 


19 


R 


13 


R 


2o 


R 


2o 


R 


2o 


4  =  cos  <#(: v  sin  a  r  ^ '  ' 

5  =   -   cos  0  (2)  cos  a  r(2) 

1  =   HlO  ,  1    R14  ,  2 

2  ~  Rj.0  (  1     Rl4  ,  1 

3  =     M.0'  3 

4  =     °14     1   f%0     5 

5  =    "*"     ^14    2    ^0    5 

R, 


-    Rg    i 


Rl5.1 


2  "   r(s) 

4  =    Rl3     1 

5  =     R-13     2 


R 


14,1 


1         "Rl°<1       r(3) 


Ris^ 

^16  ,1 

2""r(2) 

Rl8j 

4     :  :     ^14/  2    R12,4 

Rl8J 

5     =     ~     Rl4,l     Rl2    4 

^19 

1     :  :    Rl     1     R14    2 

R3,3     ^14,1 


Rl,3 

R14,l     Rl     5 
R14    2    ^1    5 
R14,] 
R14,2 
R2,4    R14.2 
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1*20,5    = 

-     R2/4    R14,l 

R2o,6     : 

R2/6 

R2L    1    : 
/ 

R3    1     R14  , 2 

R21,  2    = 

~     R3,l     R14,l 

R21,3    = 

«i.i 

Ra/  4 

R14,l     R3,5 

R21     5    : 

R14/2     R3,5 

R22,l     = 

R4#l     R14,2 

R22,2    = 

~    R4,l     R14,l 

R22    3    = 
/ 

R4    3 
/ 

R22    4    = 
i 

R14/l     R4,5 

R22,5     = 

R14,2    R475 

R23,l     = 

R5     2    R14    1 
»                  / 

R33     2    = 

1 

:    R5  f  2    R14     2 

R23,4 

:    R14     2    ^5,4 

R23,5     = 

:     "     R14,1     ^5,4 

R24    1 
f 

-         R           R}V 

■     "     R*'l      rd) 

R24,  2    : 

R24    4    _     R14,2    Rs,. 


R24    5     ~     "     R14.1     ^6,4 
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Table  2:    Transformation  Matrix  IT]  Between  Nodal  Coordinates 
and  Polynomial  Coefficients        (Columns   1-12) 


0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3T2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

T4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3T6 

0 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

3T2 

3T2 

3T2 

3T2 

0 

0 

0 

0 

0 

0 

0 

0 

3T4 

2T4 

T4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

9T6 

6T6 

3T6 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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Table  2:    (Continued)     (Columns  13  -  24) 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

Ti 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

T3 

0 

0 

0 

0 

0 

0 

0 

T5 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

2T2 

0 

0 

0 

T3 

0 

0 

0 

0 

0 

0 

T4 

0 

0 

0 

T4 

0 

0 

0 

T4 

0 

0 

0 

2T6 

0 

0 

0 

T6 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2T2 

2T2 

2T2 

2T2 

T2 

Ts 

T2 

T2 

0 

0 

0 

0 

3T4 

2T4 

T4 

0 

3T4 

2T4 

T4 

0 

3T4 

2T4 

T4 

0 

6T6 

4T6 

2T6 

0 

3T6 

2T6 

T6 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

Ti 

Ti 

Tx 

Tx 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

3T3 

2T3 

T3 

0 

0 

0 

0 

0 

3T5 

2T5 

T5 

0 

0 

0 

0 

0 
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where: 

1 


Ti   =   - 


I 


7(17 


1 

T3 =  — err 


T     = 

-•-5 


T2 

= 

1 

IrW 

= 

1 

T4 

«r" 

1 

Q£rW  i6        airW 
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IV.     NUMERICAL  CONSIDERATIONS  AND  PROGRAM  INFORMATION 

A.  GENERAL 

All  numerical  computations  were  accomplished  on  the  IBM  360/67 
digital  computer  using  FORTRAN  4  language.    The  entire  program  was 
coded  in  double  precision  arithmetic.    An  overlay  scheme  under  control 
of  the  linkage  editor  was  used  to  conserve  core  storage. 

B.  QUADRATURE 

It  was  evident  that  the  integrands  of  (3.27)  and  (3.37)  could  not  be 
manipulated  into  closed  form.    Since  the  above  expressions  were  in 
nondimensional  form,  lower  and  upper  limits  of  integration  zero  and  one 
respectively,  it  was  advantageous  to  use  Gaussian  quadrature.    An 
additional  advantage  of  using  Gaussian  forms  is  that  a  polynomial  of 
degree  2n-l  can  be  exactly  described  by  n  points  [9].    For  polynomials 
of  unknown  degree ,  maximum  accuracy  is  achieved  by  selecting  a  large 
number  of  points.    Additional  accuracy,  however,  is  obtained  at  the 
expense  of  computation  time.    Investigation  by  the  author  showed  that 
Gauss  Legendre  — 8  point  quadrature  provided  accuracy  within  the  range 
of  double  precision  round  off  error  for  shells  whose  geometry  could  be 
described  by  analytic  expressions.     For  shells  of  more  complicated 
geometry,  required  accuracy  is  obtained  when  standard  finite  element 
mesh  size  reduction  is  accomplished. 
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C .  LIMITING  CASES 

In  this  thesis,  cylindrical  shells,  truncated  cones,  and  flat  plates 
have  been  considered  as  special  cases  „    When  a  shell  is  specified  as 
one  of  the  above  forms,  geometric  approximation  is  precluded. 

Exact  solutions  for  these  shells  provided  a  powerful  tool  for  the 
analysis  of  the  stiffness  properties  of  the  general  case  of  the  curvilinear 
shell  finite  element. 

D.  CASE  WHERE  0  IS  EQUAL  TO  ZERO 

Singularities  exist  in  the  geometric  approximation  whenever  the 
element  is  constructed  in  such  a  way  that  a  point  on  its  surface  is 
locates  at  0  equal  to  zero  (see  Figure  1).     Discussion  of  the  effects  of 
this  singularity  is  contained  in  Chapter  5,  Section  B,  2. 

E.  INITIAL  CONSTRUCTION  OF  THE  STIFFNESS  MATRIX 

The  author  originally  constructed  the  stiffness  matrix  by  constructing 
all  matrices  in  (3.27)  and  then  integrating  the  entire  expression. 
Accurate  results  were  achieved.     However,  the  construction  of  the  stiff- 
ness matrix  for  a  single  element  required  500,000  bytes  of  computer  stor- 
age   and  approximately  20  minutes  of  computation  time.    This  was 
obviously  unsatisfactory  for  other  than  pedagogical  application. 

Fortunately,  identical  results  were  obtained  by  individual  stiffness 
matrix  element  computation  and  by  taking  advantage  of  symmetry.    This 
was  done  in  the  remarkably  short  time  of  15  seconds  using  only  143,000 
bytes  of  computer  storage. 
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CSFE  was  constructed  and  evaluated  in  a  separate  program  prior  to 
being  incorporated  in  a  direct  stiffness  program. 

F.      DIRECTION  STIFFNESS  PROGRAM 

Numerical  solution  of  engineering  problems  was  accomplished 
through  the  use  of  a  stress  analysis  program  which  was  based  on  the 
direct  stiffness  method.    The  total  stiffness  matrix  and  the  load  vectors 
were  assembled  using  programming  techniques  furnished  by  Professor  G. 
Cantin.     This  program  gave  numerical  values  for  displacements ,  moments, 
stresses,  and  strains  at  each  nodal  point  of  the  mesh.    A  flow  diagram  of 
the  computer  program  is  presented  in  Figure  3.    It  can  be  seen  that  the 
program  is  functionally  divided  into  two  parts.    In  the  CSFE  portion,  the 
stiffness  matrix  and  its  associated  pressure  load,  strain-displacement, 
and  stress-strain  relationships  are  computed  for  individual  elements. 
Results  of  these  computations  are  stored  on  disk  units.    In  the  stress 
analysis  portion,  individual  element  information  is  retrieved  from  disk, 
and  assembled  in  a  manner  consistent  with  the  direct  stiffness  method. 

The  author's  primary  purpose  for  construction  of  the  stress  analysis 
portion  of  the  program  was  to  verify  the  stiffness  properties  of  the 
curvilinear  shell  finite  element.    It  was  therefore,  designed  to  handle 
problems  of  limited  size.     Specifically,  the  program  is  limited  to  a 
maximum  of  36  elements,  49  nodal  points,  or  294  equations.    Coding 
methods  used,  however,  are  general  to  allow  for  future  program  size 
expansion.    The  computer  program  is  found  in  Appendix  C. 
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V.     DISCUSSION 

A.  INTEGRITY  OF  THE  CURVILINEAR  SHELL  FINITE  ELEMENT 
Confidence  in  the  integrity  of  CSFE  was  established  by  comparing 

stiffness  properties  of    limiting       cases  (flat  plate  and  cylindrical  shell) 
with  properties  found  in  the  literature.    Element  by  element  comparison 
with  values  given  by  Bogner,  et.  al .   [Ref.  10]  revealed  exact  agreement 
for  the  case  of  the  flate  plate.    In  the  case  of  the  cylindrical  shell 
element  similar  comparison  was  made  with  values  given  by  Cantin  [Ref. 
1];  once  again  exact  agreement  existed. 

Confidence  in  the  integrity  of  the  stress  analysis  program  was 
established  by  direct  comparison  of  problem  solutions  using  CSFE  with 
solutions  found  in  the  literature.    Results  of  comparison  are  discussed 
in  the  following  sections. 

B.  NUMERICAL  SOLUTION  OF  CLASSICAL  PROBLEMS 

The  classical  problems  analyzed  in  this  section  were  selected  to 
demonstrate  the  integrity  of  the  element,  and  to  exhibit  both  the 
flexibility  and  limitations  associated  with  CSFE. 

1  .      Pinched  Cylinder 

The  pinched  cylinder  illustrated  in  Figure  4  has  been  analyzed 

by  Bogner,  et  al    [Ref.   11]  using  a  (48x48)  element  and  by  Cantin 

[Ref.  6]  using  a  (24x24)  element.    Table  3  is  a  comparison  of  results 

obtained  using  CSFE  with  those  found  in  the  literature.    The  obvious 

agreement  of  results  confirms  the  integrity  of  CSFE.    Moreover,  it  is 
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Fig.  4  PINCHED  CYLINDER 
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evident  that  CSFE  has  the  desirable  characteristic  of  rapid  convergence 
to  the  true  solution  (-0.1139)  given  by  Cantin. 


No.  of        Bogner    et    al  Cantin  CSFE 

EQS.  (48x48)  (28x28)  (24x24) 

Mesh           Displace.       Mesh     Displace.     Mesh      Displace. 
(inj 

54  2x2  -0.0931         2x2  -0.0932 

108  2x2  -0.0808 

150  4x4  -0.1126         4x4  -    .1128 

180  2x4  -0.1098 

294  6x6  -0.1137         6x6  -0.1138 

486  8x8  -0.1139 

726  10x10  -0.1139 

Table  3.    Displacement  under  load  P  for  pinched  cylinder 

2  .      Spherical  Cap 

The  spherical  cap  illustrated  in  Figure  5  demonstrates  the 
general  case  of  CSFE.    This  problem  which  has  an  exact  solution  given 
by  Timoshenko  [Ref.  12]  is  useful  for  demonstrating  the  capabilities  and 
limitations  of  the  new  element.    Figure  6  contains  a  plot  of  meridional 
bending  moments.2 


2  The  signs  of  Timoshenko' s  results  have  been  changed  to  agree  with 
the  author's  convention. 
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It  is  worthy  of  note  that  the  pole  of  the  cap  was  located  at  a  singular 
point  (0  =  0).    In  order  to  fulfill  connectivity  requirements,  a  small  rigid 
polar  "cap"  was  assumed  to  exist.    Boundary  conditions  were  applied  to 
allow  this  "cap"  freedom  of  displacement  in  the  z-direction  at  zero  slope. 
Zero  slope  specification  had  the  effect  of  introducing  a  false  moment  in 
elements  adjacent  to  the  polar  cap.    This  false  moment  is  visible  in 
Figure  6.    It  is  important  to  recognize,  however,  that  the  detrimental 
effect  of  the  erroneous  specification  is  localized  in  the  near-pole 
elements . 

Results  indicate  that  the  bending  predictions  are  quite  accurate 
even  though  a  coarse  mesh  was  used. 

C.     CONCLUSIONS 

A  finite  element  for  a  shell  of  revolution  which  has  arbitrarily  curved 
meridians  has  been  developed  and  tested.    Element  description  was  based 
on  displacement  functions  which  met  all  requirements  to  guarantee  mono- 
tonic  convergence  to  an  exact  solution  by  mesh  size  reduction. 

Results  for  classical  problems  obtained  by  using  the  new  element 
indicated  that  good  approximations  for  bending  stresses  could  be  obtained 
with  a  very  coarse  mesh.    Mesh  size  reduction  will  be  necessary  to  judge 
element  performance  in  problems  where  membrane  action  is  predominant. 

The  geometrical  deficiency  where  ($  =   0)  is  a  limitation  which  must 
be  recognized  by  users.    Its  effect  is  problem  dependent.    For  the  case 
of  the  spherical  cap  it  was  seen  that  the  detrimental  effect  of  artificially 
imposed  boundary  conditions  at  this  singular  point  did  not  propagate 
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throughout  the  solution.     On  the  other  hand,  the  author  was  unable  to 
successfully  approximate  stresses  in  a  pressurized  toroidal  shell  due  to 
error  propagated  from  similar  artificial  boundary  conditions  .    Time 
precluded  the  determination  of  a  method  to  efficiently  handle  this 
singularity. 

Not  withstanding  this  geometric  limitation,  results  indicate  that  the 
new  element  (1)  is  accurate,     (2)  requires  low  computational  effort  and 
(3)  does  not  require  an  analytical  description  of  shell  geometry.    Based 
on  these  facts  it  is  concluded  that  the  curvilinear     finite  element  can  be 
efficiently  used  in  many  engineering  applications. 

D.     RECOMMENDATIONS 
It  is  recommended  that: 

(1)  the  stress  analysis  program  be  expanded  in  order  to  handle 
realistic  engineering  problems; 

(2)  a  consistent  gravity  load  vector  be  included  in  the  stress 
analysis  program;  and 

(3)  a  study  be  conducted  to  determine  a  method  of  efficiently 
handling  the  singular  point  which  exists  where  (0  =   0)   in  the  new  element 
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APPENDIX  A    NUMERICAL  EXAMPLE  OF  THE  STIFFNESS  MATRIX 
In  this  appendix,  a  typical  example  of  the  stiffness  matrix  CSFE 
is  presented.    Eigenvalues,  before  and  after  rigid  body  motion  was 
included,  and  matrix  K22  are  also  presented.    Double  precision  figures 
were  truncated  to  conserve  space. 
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STIFFNESS    MATRIX    FIR     t    TYDICH     SHf-LL    ELr.MENT 


YCljNG'S    MODULUS    =  1  .OOOOOO 

PCISSnMS  WiTIT    =  0.300000 

RADIUC    1  =  7.071100 

RADIUS    2  =  8„,660100 

L^NOTH  =  2.MP.000 

DHI     1  =  0.7^5^00 

PHI    2  =  1.0^7200 

TH!C*N=SS  1                =  1.000000 

THICKNCSS  ?               =  1. 000000 

ALPHA  =  0.?M*00 
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■02-1.38200- 
•03-1.72910- 
■03-7.13170- 


•03  1.40020-02 
•02  1  .59000-02 
•01  5. 019^0-02 
•02-2.50830-02- 

01  4.357*0-02 

02  1.26150-01 
02-2.51250-02- 
02    1.673  c0-02 
02    3.1°620-02 

0  2    *-.1215  0-0  2 
0?    3.39070-02 
02-3. 03520-02- 
O2-2.39S50-O2- 
02-1  .2751  C-02- 
02    3.75500-02- 
02    2.639rD-0?- 
02    7.05530-0"* 
03-5.55010-0  J- 
0->    2.^4990-02 
02-2.32040-0  2- 
01-1  .33200-01- 

01  1.52680-01 

01  8.89860-02 

02  6.10^70-02 


3.43610- 
3.  V4840- 
1.31320- 
4. 357*0- 
?.  09810- 

2.°5qon- 

3.3  ?nct0- 
2  .  7  6  8  4  T  - 
3.108br- 
3.'  3C30- 
4  .70240- 
3.1572r- 
^.a112C- 
5.^6OS0- 
1.01  37C- 
1.06240- 
2.03770- 
1.88300- 
4.0493C- 
6. 53560- 
1.72910- 

1.73410- 

5.^5350- 


•02  1.54690-02 
03  1.28580-03 
•01  3.C674C-03 
02-1.26150-02 
02    2.959C0-02 

02  5.84?2P-0-> 
02-1  .65570-02 

03  2.56Q40-0'' 
02  2.8087^-02 
0?  2.40430-02 
02  2. 76*10-02 
02-1  .68880-02 
02-1. 73170-02 
04-1 .97280-0^ 
02-5  .401«r-0? 
03-3.31520-0^ 
02  2.10610-02 
02-1  .80960-02 
02  1.78440-02 
03-1.9920r-03 
01-7.1 3170-02 
02  6.10*70-0? 

01  6.46350-02 

02  4.6593C-02 


(COL.  lc  7C  2^  ) 
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£IGr-NVALLcS    RPFTRF    PIGIC    BCDY    "CTICN    INCLUDED 

1.6C7<;r    CC    8.93580-01    8.03190-01    7.5090D-01    S. 32120-01  c. 91280-01 

4.474CD-C1    4.12140-01    1.68540-01    1.01600-01    9.7839D-02  8.4<?P60-0? 

f#1Cccn-C2    4.83840-02    1.87230-02    1 .  IS40* H-O7    1.72800-02  f. 22710-03 

4.516CD-C2    ^.e^Q^r-O^1    2.62960-04    3.ci">45n-05    1.09450-06  1.5^470-07 


VATRT *  K22 

2.15120-C4-i.64  9cn-OR-l  .64980-0*  6  .  =J63 -S  0-0  c  1  .69770-04- 2.  0<»39C-05 
4. **<;<;  D-Cc  c;.621,;[:-04-2.17210-0=;-6.5">&c>D-0'f-6.i;6  36C-05  1  .«;905C-04 
1  .6<.TcO-C^-:.1721D-C5  2.14*40-04-1  .036  20-04  8  .  2^0^  0-0*-o  .  605  20-  1  7 
5.56  36D-C5-6.5969C-04-!  .08*2  0-04  5  .  3  24  3  0-0  3-3  .^otl  3D-0<>-6  .  16  240-0  3 
1.69770-04-6.56360-05  8.  25  0  80-04-^  .6^3  3  0-0 '♦  1  .2  3  17r_0?-3  .  1 1  30C-0* 
2.09  390-C5  1.59C5r-0^-2.63630-]  5-6  .  1624-^-0  3-8  .  1 1  30  0-04  8  .  cO"?*o-o^ 


FIGENVALtFS  ftPTFR  RIGID  PQDY  MOTION  WAS  TNCLUCEO 

1.6C3CD  CC  8.7532C-01  7. 03650-01  5.7*820-01  =1.70550-01  4.6P02D-01 
3.765er-Cl  2.86UD-01  1.03530-01  1.00830-01  3.->773r-02  7.2P31D-0? 
4.2277D-C2  ?. 06*90-02  1.7525D-02  1.50520-02  L.1759C-02  c. 32370-03 
9.2^71C-!7  C.1124D-17  3.91360-17-7.79230-18-4.65150-17-7.73010-17 
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APPENDIX  B    COMPUTER  OUTPUT 
In  this  appendix ,  a  complete  set  of  results  for  the  spherical  cap 
analyzed  in  this  study  are  presented.    Results  for  the  spherical  cap  were 
obtained  using  a  4x4  mesh.    The  mesh  layout  accompanies  the  problem. 
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SPHFRICAL  SHELL  WITH  EXTERNAL  PRESSURE  (1  PSI) 

TOTAL  NUMBER  ^p  ELEMENTS 

M1MPFR  OF  0IFFFo»=NT  TYPE  ELEMENTS 

TOTAL  NUMPER  CF  NODES 

NUMBER  nc  NODES  W!TH  BOUNDARY  CONDITIONS 

HALFRAN'D  UID^M 

NUMBER  Or  POINTS  WITH  CONCENTRATED  LOADS   = 


16 

2C 
2* 

42 

0 


I  6  I.I  16  21 


MESH  LAYOUT 
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FL 

SPEC 

f  /nris 

RA01/PAP? 

PHII/PHI2 

1 

REN 

A. 00000  06 

^.1^10  00 

3.4910D-02 

1  .66700-01 

1.56280  01 

1.74530-01 

2 

r,FN 

4.00000  Of. 

1.  .56  2  30  0! 

1 ,71Kin.0' 

1  .66700-01 

?.3C570  01 

3.22990-01 

HI  /H2 


L     /HPHA 


3.0000D    00       3.~27C0-C1 

3.noccn   cr     i. 13510  ri 

3.00000    00       3.!;27Cr-Cl 


GEN  ^-.COCOO    06       2.«R57n    01       3. '2890-01        ^.OOOCC    OC       1.3351p    0! 

1.66"'00-01       A.OR5oo    o1       4  .7 1  24D-0 1       ^.OODCn    00       3.927CI>C1 


^PN  '.00000    06       '.08510    01       '.712<.0-0l 

i. 667QD-0!        C.1622P    0!        6.10870-01 


*  .00000  00   1  .2c  +  7~    0! 

^.oooon  oo  3.«?2"'cr-r1 


CONNECTIVITY 

T  J 

1  2 

?  3 

^>  4 

'  5 

f  7 

7  3 

P  0 

c  10 

11  12 

12  13 
1  3  14 
14  15 

17  19 

IF  19 

ic  20 


7 

8 
1 

'0 
12 
13 
14 

15 

17 
!  8 
i  o 

20 
22 
23 
24 


a 

Q 

i : 

12 

13 
14 
i  « 
i  -i 
1  B 

1  = 
21 
22 

2  3 

76 


Type 
i 

■> 
■x 
.'- 

1 

2 
-j 
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RRCNOAPY  CONDITIONS 


JOINT 

X     TRANS 

Y    TRANS 

I    TRANS 

w,  x 

W,Y 

w.xv 

1 

0 

0 

0 

0 

c 

f~ 

0 

0 

n 

0 

c 

11 

0 

0 

0 

0 

c 

16 

0 

0 

0 

0 

c 

21 

0 

0 

0 

0 

c 

5 

0 

0 

0 

0 

0 

c 

1C 

0 

0 

0 

n 

0 

c 

15 

0 

0 

0 

0 

0 

c 

20 

0 

0 

0 

0 

c 

c 

2  = 

0 

0 

0 

0 

0 

c 

2 

0 

1 

0 

c 

i 

0 

1 

0 

c 

u 

0 

1 

0 

c 

■>? 

0 

1 

0 

c 

2? 

0 

1 

0 

c 

24 

0 

I 

0 

c 

7 

1 

c 

0 

P 

1 

0 

c 

c 

1 

0 

c 

12 

1 

c 

c 

13 

1 

c 

c 

14 

\ 

0 

c 

17 

1 

0 

c 

1? 

1 

n 

c 

!=• 

1 

0 

r 

0!  SPL 

?r"*cMS 

JOINT  L 

1  -c.c 

2  ?. 09120-05 

3  f.4?3CD-0'5 

4  fc.i«?4D-CB 

5  -C.C 
ft  -CO 

7  3.  C5120-05 

Q  £. 1*240-05 

10  -c.c 

11  -c.c 

12  3.0312D-CK 

13  6.4P3C0-05 
1*  6.1=240-0= 
1  =  -C.C 

1*  -C.C 

17  3.CH2O-05 

IP  6.4?3C0-0= 

1=  6.  1=2  40-0*; 

2C  -C.C 

21  -C.C 

22  ?. 0^120-0  = 

23  6.*83CD-0= 

24  *-. 1  =  240-0  = 
?  =  -C.C 


V 

W 

W,X 

W,Y 

O.C 

-6.00219-0* 

-0.0 

-o.c 

-0.0 

o.c 

-5.6M10-04 

-6.i?5o^-06 

-c.c 

-0.0 

o.n 

-i.rjp-o' 

-I  .50350-05 

-0.0 

-0.0 

0.0 

-1.R2260-04 

-?  .15830-05 

-0.0 

-o.c 

0.0 

-0.0 

-0.0 

-0.0 

-0.0 

o.c 

-6. 00210-0' 

-0.0 

-o.c 

-co 

3.20510- 

-17 

-K.f«.i  10-04 

-5  .12  =  00-06 

-0.0 

-co 

?  .C4<»"0- 

•1« 

-^.323n-04 

-1  .50350-0=. 

-O.o 

-co 

l  .759=0- 

■1ft 

-1.82260-0' 

-2.15"3D-05 

-0.0 

-o.c 

0.0 

-0.0 

-0.0 

-0.0 

-o.c 

0.0 

-6.00210-04 

-0.0 

-0.0 

-0.0 

4.8160D- 

■1T 

-5.6-1 1^-04 

-5  .12500-06 

-0.0 

-o.c 

<■  .30830- 

-17 

-4. 32310-04 

-1  .50350-0  = 

-co 

-c.c 

2.2^*20- 

■IS 

-1.82260-04 

-2.1=P1O-05 

-0.0 

-o.c 

0.0 

-0.0 

-0.0 

-co 

-c.c 

0.0 

-6.00210-04 

-0.0 

-co 

-o.c 

3.20740- 

-17 

-5.64i in-04 

-5.12500-06 

-co 

-o.c 

1  .99380- 

-11 

-4.32310-04 

-1  .50350-0  = 

-0.0 

-0.0 

1.7593D- 

-IS 

-1 .8226n-o<- 

-2.15R3D-05 

-0.0 

-0.0 

O.C 

-0.0 

-0.0 

-co 

-co 

0.0 

-6.00210-0* 

-0.0 

-co 

-o.c 

0.0 

-5.6M  10-04 

-5.12500-06 

-0.0 

-0.0 

0.0 

-i. 32310-04 

-1  .503  =  0-05 

-co 

-co 

0.0 

-1.82260-04 

-2. 15830-05 

-co 

-cc 

0.0 

-0.0 

-0.0 

-co 

0.0 

W,XY 
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STRF5CM 


JCTNT       NX 

i  -1.12990  02 

2  -4.s4i7o  01 

3  -5.0434P  01 

4  -=.p'3co  oi 

=  -6.C429D  01 

*  -l'.  32990  02 

7  -4.04170  01 

P  -=.C434P  01 

0  -5.P23CD  01 

10  -6.C429D  01 

11  -1  .  3705D  02 

12  -4.94170  01 
1?  -=.C434D  01 
1^  -=.°23CO  01 
1  "S  -t,c»?n  oi 
15  -1*32990  02 
17  -4.04170  01 
1"  -=.043*0  CI 

lc  -6.  °2,cn  01 

20  -6.04250  oi 

21  -l."3?eeP  02 
2?  -4ie*l"»D  01 

23  -=.04340  0! 

?.A  -C.<>23C0  01 

15  -6.10*290  01 

JCIN1      SNX 

1  -r.66rFD-C6 

2  -3.2936P-C5 
»  -3. '4420-05 
'  -4.6C370-06 

=  -4.e<;^?o-c6 

'  -C.665SD-C- 

7  -?.283tn-G6 

0  -3.64420-06 

0  -4. 60^70-06 

]C  -4.8CCO0-05 

il  -c. '66eo-06 

12  -•>.2Q,60-06 

13  -3.64420-06 
'.4  -4.60370-06 
15  -4.P9=P0-06 
I*  -9.665PO-06 
l7  -3.23360-0'. 
IP  -3.64420-05 
1Q  -4.6C370-05 

20  -4.P0500-06 

21  -c.t66ep-06 

22  -3.29360-0* 
2">  -3.64420-06 
74  -4.fcC370-06 
?e  -4.PO  =  fD-0iS 


NY 

-1.0200D  02 
-6..  00750  01 
-4.02160  01 
-1.79080  01 
-1.00730  01 
-1 .02000  02 
-t>. 00750  01 
-4.02160  01 
-1.7O080  01 
-1.00730  01 
-1.02000  0? 
-6.'0075D  01 
-4^02160  01 
-1 .7=030  01 
-1.00730  01 
-1.02000  02 
-6.C075D  01 
-4.02160  01 
-1.79080  01 
-1.0073O  01 
-1.02000  0' 
-6.00750  01 
-4. 021*0  01 
-1  .7°080  01 
-1.00710  0! 


NXY 

0.0 
-3.80500-12 
-3.86360-1 

-5.10700-1 

0.0 
-1.31160-1 
-1.3746D-1 

7.35530-1 
-4.47280-1 
-7.2004O-1 
-1.97090-1 
-2.11660-1 

2.00610-1 
-3.07170-1 
-0.30660-1 
-1.31250-1 
-1.37630-1 

7.32590-1 

-4.474'0-l 
-7.1O970-1 

0.0 

-3.81210-1 
-3.87030-1 
-5.111=0-1 

0.0 


"X 

-4.68020  00 
-4.44680  00 
-7.62730  00 
-2.61650  00 

3.  17O00  01 
-4.5«o?0  00 
-4.  .46^0  00 
-7.62730  00 
-2.61650  00 

3.17QSD  01 
-4. 68=20  00 
-4. .4680  00 
-7.42730  00 
-2.61650  00 

2  .  17900.  01 
-*.*3920  CO 
-4.44*30  00 
-7.62710  00 
-2.61650  00 

3.17QO0  oi 
-4.6=020  00 
-4.44680  00 
-7.62730  00 
-2.616=0  00 

3.  17OQ0  01 


MY 

-7.8170D-01 
-3.46810  00 
-5.55020  00 
-4.53750  00 

=.30100  00 
-7.P1700-01 
-3.46M0  00 
-c. 55020  00 
-4. 5^7  =  0  OC 

=.^0100  00 
-7.P17C0-01 
-3.46810  00 
-=.==020  00 
-4.517=0  00 

5.301  CO  00 
-7,ei7CD-0] 
-3.45810  00 
-5.55020  00 
■*.537=o  00 

c. 30100  00 
•7.P170O-01 
■i.  45810  00 
■5.5=020  00 
■'  .  =  3750  00 

=  .301C0  00 


"XV 

0.0 
-2.02540-14 
-3.2112r-' 3 
-4.27010-1 3 

0.0 
-3. 12380-12 

8.1  1*40-1  3 

■7 .4,4Bir-i  1 

-9.50480-12 
-1.C461C-12 
-4  .601CO-! 2 

1  .62'20-!2 
0  .Cl 770-1? 

-1  .C  =  740-1  2 

-I ,36=?r-l  2 

-3.12610-1 2 

a. 11 °C0-1  2 

7.645*0-1 3 

-  =  .r  °=cr-: 3 

-1 .C45C0-1 2 

0.0 

-2.01] 'P-l* 

-3.21*60-13 
-J4. 27*30-1? 

1 .200=0-2  = 


SNY 

-6. 6=250-06 
-4.  3,1980-06 
-2.  '=070-06 
-6.834*0-07 
2.90520-22 
-6.6=250-06 
-4.31780-06 
-2.65070-06 
-6.8344D-07 
-1.3  280O  -2' 
-6.65250-06 
-4. 31O8O-06 
-2.6507.0-06 
-6  .  93440-07 
-3.32500-22 
-6  .56?5D-06 
-4.31980-06 
-2.65070-06 
-6.93440-07 
-3.32«00-->2 
-6.6=250-06 
-4.31 980-06 
-2.6607D-06 
-6.83440-07 
-9.6=130-27 


SNXY 
0.0 
-7. 30670-10 

-"•.51320-1P 
-9.93040-18 

0.0 
-2.5=0=0-18 
-2.67300-18 

1  .43010-18 
-3.6°740-lP 
-1  .40010-17 
-3.8324D-1P 
-4.11500-1° 

3.90040-1= 
-5.07110-10 

-1  .82720-17 
-2.5=210-18 
-2.67670-18 

1.42300-1P. 
-3.70110-1P 
-1  .40000-17 

0.0 
-T.4l 370-1C 
-7.52610-18 
-9.93990-18 

0.0 


S*X 

-=  .06  =  =  0-07 
-4.23P60-07 
-7. \4(  70-07 
-2.06680-07 

3. -.3510-06 
-5  .05550-07 
-4.7JP50-07 
-7.44670-07 
-2.06680-07 

3.41510-0* 
-5.06550-07 
-4.2 5  0=1-07 
-7. 44670-07 
-2.06680-07 

3.-t351O-06 
■5  .06*60-07 
-4.298=0-07 
-7.44670-07 
-2.06630-07 

3.43510-0' 
-5.0o6rP-07 
-4.  708  =  0-07 
-7. '.4670-07 
-2  .06^-80-07 

3.43810-06 


SMY 
-2.' 6840-1= 
-?.  02900-0"" 
-'-."'6410-07 
-4  .  =  =  700-07 
-3.^410-2' 
-1  .10230-1P 
-1. 02080-07 
-4. "'5410-07 
-4.==7CD-07 
-=.03680-22 
-1  .  =0740-19 
-3  .C'980-07 
-4.  -"5410-07 
-«.  =  =  700-07 
-3.036PO-22 
1.27300-1 3 
-3. 029OO-07 
-4 .75410-07 

-* .6570O-07 

-1. 03680-22 
3.70470-1 9 
-3.02OPP-0"7 
-4.75410-0"' 
-t  .  =  =  700-0'' 
-2.32960-22 


SMXY 


O.C 


-7.5.B*4n-21 
.p 1  ^^  2  =  40- 20 
-1.1 C72T-1 0 

O.C 

-  3 . Cc-0O-i c 

2.1 C«30-ia 

l.9°3C0-19 

-2.4 "tpp-] o 

-'."122^-10 

- 1 . 1 1  7C0-1  8 

-»  .  3CAcr_i  : 

' . C7 " °P_ ; c 

-'.s*c20— 1 c 

. 1 . 5 1 c  *n- 1 q 
-0. 1C500-1 o 

2.1C470-1C 

I  .ca-Jin-'  o 
-2.4  8=10-' Q 
-2.71 20D-1 0 

0.0 

-7.*CC3n-21 
i9. 3*960-20 
-I  .1 C?r-i c 

3.i'"2r-i2 
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APPENDIX  C    COMPUTER  PROGRAM 
In  this  appendix  the  stress-strain  analysis  program  is  presented. 

User  and  functional  information  are  presented  through  the  use  of  comment 

cards  located  throughout  the  program. 

The  program  input  data  deck  is  illustrated  below.    Data  values 

presented  were  those  used  for  the  spherical  cap  (2x2  mesh). 
12  345  6789012345678901234567  89012345678901234567890 


'/* 

Orange 

card 

- 

rr~ 

BLANK 

r~ 

BLANK 

/ 

BLANK 

f 

9 

f 

8 

1 

1 

/ 

7 

6 

<r 

5 

1 

1                     1 

1 

r 

\ 

4 

/ 

3 

/ 

2 

1 

1 

/ 

1 

/ 

5 

6 

9                    8 

2 

/ 

4 

5 

8                    7 

1 

2 

3 

6                    5 

2 

/ 

1 

2 

5                    4 

1 

^00000. 

.1667 

3. 

3.                  -1. 

'2  7 

.8118 

51. 

6222 

27. 

6039         .31416          .61087 

.7854 

'4000000. 

.1667 

3. 

3.                  -1. 

1 

1491 

27. 

8118 

25. 

,1325         .03491         .31416 

.7854 

/ 
i 

4- 

TEST 

4 
PROB 

2 

.  (spheric; 

9                    9 

30 

\L  CAP    2x2  MESH) 
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C   STRESS  ANALYSIS  PROGRAM  FOR  CONSTRUCTING  AND  TESTING  CSFE 

C 

C   CONTROL  SUBPROGRAM 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON/Bl/c,Hl,H2,POIS,RA01,RAP2,L,PHI  1 , PHI?, ALPHA, 
#THETA,KASE 
COM VON/MB N/8 LOAD (294) , NBC (50,7) ,NCON(36,5) 
C OMMON/P ARM /NEL, NELT, NDPT, NPBC, NBA ND,NFQ 
50  DO  100  1=1,294 
100  BLOAO( I )=C.OOO 
CALL  INPUT 
IF(NELT.EC.O)  STOP 
CALL  OV1 
CALL  0\J2 
GO  TO  ^0 
ENn 

SUBROUTINE  0V1 
IMPLICIT  REAL*R  (A-H,0-Z) 

COMMON /MRN/BLOA 0(294) ,NBC( ^0  ,  7  )  ,  NCON  (  36  ,  5) 
C0MM0N7PAPM/NEL , NELT , NDPT , NPBC , NBAND, NFQ 
DIMENSION  TDTK( 294,54) 
DO  100  1=1,294 
DO  100  J=l,54 
100  TOTK(  I  ,J)=C.CDO 
CALL  FORML 
CALL  MERGER  (TOTK) 

CALL  BANSOL( TOT K,B LOAD, NEO,NB AND) 
RETURN 
END 

SUBROUTINE  0V2 

IMPLICIT  REAL*8  (A-H,0-Z) 

CALL  CUTPUT 

CALL  PSSN 

RETURN 

FNO 

SUBROUTINE  INPUT 

IMPLICIT  RFAL*8  (A-H,0-Z) 

REAL*8  K1,K2,K?,L 

INTEGFP*4  SPTY(4),SPEC 

DIMENSION  STFM0D(24,24) , T I NV ( 24 , 24 ) , TI NVT ( 24, 24) 

DIMENSION  f>VEC(24)  ,CF0RCE(6) 

C0MM0N/P1/E,H1,H2,P0IS,RAD1,PAD2,L,PHI 1 , PH I  2  ,  AL PHA , 
#THETA,KASE 

C0MM0N/MBN/BL0AD(294)  ,NBC(50,7)  ,NC0N(36,5) 

COM MHN/P ARM /NEL, NELT, NDPT, NPBC, NBA NO, NFQ 

COMMON/SS  SN/SS(24, 24) ,SN( 24,24) 

DATA  SPTY/»GEN  CYL  PLT  CONE'/ 
C   BLOAD  IS  TOTAL  LOAD  VECTOR 

C   NBC  IDENTIFIES  THE  NODE  POINT  WITH  ITS  BOUNDARY  CONDITION 
C   NCON  IDENTIFIES  ELEMENT  NODAL  POINT  AND  ELEMENT  TYPE 
C   NBAND  IS  ECUAL  TO  HALFBAND  WIDTH 
C   NPCL  IS  THE  NUMBER  OF  POINTS  WITH  CONCENTRATED  LOADS 

NTK  =  0 
C   NTK  IS  CYLINDER  TRACK  COUNTER 

WRITE(6,1100) 
1100  FORMAT( 1H1,//) 

READ  (5,1000) 

WRITE  (6,1000) 
1000  F0RMAT(80H 

n  ) 

READ  (5,2000)  NELT , NEL , NDPT , NPBC , NBAND, NPCL 
2000  FORMAT(MIO) 

IF(NFLT.EC.O)  RETURN 
NEQ=6*NDPT 
C   NELT  IS  TOTAL  NUMBER  OF  ELEMENTS 
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NEL  IS  NUMBER  OF  DIFFERENT  TY 
NDPT  IS  TOTAL  NUMBFR  OF  NODES 
NPRC  IS  NUMBER  OF  NODE  POINTS 
WRITF  (6,2500)  NELT,NEL,ND 
2500  FORMAT* /,39X, • TOTAL  NUMBER 
#38X,  'NUMBE*}  CF  DIFFERENT  T 
#38X, 'TOTAL  NUMBER  OF  NODFS 
*38X, 'NUMBER  OF  NODES  WITH 
#/,38X,  'HALPBANn  WIDTH', 28X 
JMNUPBER  CF  PCINTS  WITH  CON 
WRITF  (6,3000) 
^000  FORMAT (1H!  ,///,32X,* EL    S 
**RA02    PHI1/PHI2      Hl/H 
DO  100  1=1, NFL 
READ(5,^500)  R AD1 , RAD2 , L , P 
F0RMAT(7F1C. 0) 
READ  (5,4000)  F,P0ISfHl,H2 
F0RMAT(5F10. 0,110) 
IS  THE  PRCSSUPE  LOAD 

SPECIFICATIONS  FOLLOW 

KASF=0 

KASE=1 

KASE=2 

KASF=3 

S  SPECIFIED 

ALPHA=1.0D0/ 

SPEC=SPTY(1) 

SPEC=SPTY(2) 

SPFC=SPTY(3) 

SPEC=SPTY(4) 

T,SPEC, E,RA 


PE    ELEMENTS 


WITH    BOUNDAPY    CONDITIONS 
PT,NPBC,NBAND,NPCL 

OF  ELEMENTS*  ,18X,* =•  ,  110,/, 
YPE  ELEMENTS* ,PX,*=« ,110,/, 
K9iv.isi.fin./. 


YPE  rucci' 

•  ,nx,*  =  * ,  no,/, 

BOUNDARY  CONDITIONS 
, *=*,  I10,/,38X, 
CENTPATED  LOADS 


=  '  »no, 

=', 110) 
RAO!/', 


3500 


4000 


P 
KASE 


IF 


A  FLAT  PLATE  I 
IF(KASE.EQ.2) 
IF(KASF.EO.O) 
IF(KASE,cQ.l  ) 
IF(KASr.FQ.2 ) 
IFCKASE.EQ. ^) 
WRITE  (6,4100) 
#H2, ALPHA 
4100  FORMAT(//,32X,I ?,3X,A4,2X, 
CALL  ELSTF(STFMOD) 
NTK=NTK+1 

CALL  WRDI SK(NTK,STFM0D,576) 
CALL  TFCRW(TINV,TINVT) 
CALL  STRESS(TINV,SS,SN) 
NTK=NTK+1 

CALL  WRCISK(NTK,SN,576) 
NTK=NTK+i 

CALL  WRrTSK(NTK,SS,576) 
NTK=NTK+1 

CALL  PLCAD(P,PVFC) 
CALL  WRriSK(NTK,PVEC ,24) 
CONTINUF 
WRITE  (6,4400) 


PEC      E  /POIS 
?       L  /ALPHA*  ) 

HI1,PHI2, ALPHA 


CURVILINEAR  SHELL 

CYLINDRICAL  SHELL 

FLAT  PLATF 

TRUNCATED  CONE 

WIDTH  IS  WRITTFN  FOR  ALPHA 

THETA 


D1,PHI1,H1,L,PPIS,RAD2,PHI2, 
IP 50 12. 4// 43 X, IP 5D 12. 4) 


100 
4400 
450C 


FORMAT ( 1 HI ,///f32X, 'CONNEC 

WRITF  (6,4500) 

FORMAT ( //,33X,*I*,15X,*J*, 

DO  200  I=1,NELT 

READ  (5,5000)  (NCON(I,J),J 

FORMAT( 51 10) 

WRITF(6,5500)  (NCON(I,J),J 

FORMAT ( /, 32 X, 12,14 X, 12,1 4X 

CONTINUF 

WRITF  (6,6000) 

FORMAT ( 1  HI ,///, 2 7X,* BOUNDARY  CONDITIONS*  ) 

WRITE(6,6100) 

F0RMAT(//,27X,*J0INT',3X,» 


SOOO 

5500 
200 

60C0 

61C0 

#*Z  TRANS1 ,8X,*W,X* ,9X, • 
READ  BOUNCARY  CONDITIONS 
DO  300  I=1,NPBC 
READ  (5,65C0)  (NBC(I,J) 
6500  F0RMAT(7I10) 

WRITF(6,7000) 
7000  FORMAT!/, 29X, I 
#11X,I1) 
300  CONTINUE 

IF(NPCL.EQ.O) 
WRITE  (6,7050) 
7050  F0RMATUH1,///, 


TIVITY' ) 

15X,*K*,15X,*L*,13X,*TYPE*) 

=  1,5) 

=  1,5) 
,I2,14X,I?,1^X,I?) 


X  TRANS*, 5X,*Y  TRANS1,  5X, 
W,Y*  ,9X,  *W,XY*  ) 


J=l,7) 


(NBC(I,J) ,J= 
2,  8X,I1,UX 


RETURN 


1*71 

,11 ,11X,I1  ,11X,I1  ,11X,I1, 


27X, 'CONCENTRATED  LOADS*) 
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•3X,»F0RCF  XSSXt'FDPCF  Y»,5X, 
•  W,Y'  ,<?X,  'W,XY«) 


WRITE  (6,7100) 
7100  FORMAT (//,27X, • JOINT 
#»FORCE  7*  ,RX,«W,X»  fPi 
READ  CONCENTRATED  LOADS 
DO  500  I=1»NPCL 

READ  (?,72C0)  JNR  ,  (C.FORCE  (  K  )  ,K  =  1  ,6  ) 
7700  FORMAT  (I10,6F10.0) 

WRITF  (f,7300)  JNR, (CFORCE(M) ,M  =  1,6>) 
7300  F0RMAT(/,19X,I 12, 1P6D12.4) 
IR=6*JNR-6 
DO  400  J=l,6 
IRL=IB+J 
400  BLOAD( IPL)=CFOPCE( J) 
500  CONTINUE 
RETURN 
END 


FUNCTION  COSPHI (X) 
IMPLICIT  RE4L*8  (A-H,0-Z) 
REAL*0  K1,K2,K3,L 

COMMON /R 1 /F, HI, H2 f PO TS,R AD  1, RAD? tLt PHI  1 , PH I  2 , ALPHA, 
*THETA,KASE 
C0MM0N/R2/CPHI1 ,CPHI2,C0,C1,C2 
COSPHI      =  C0+C1*X+C2*X**2 
IFIKASE.GE.l )   COSPHI=0.000 
IF(KASE.EQ.3)   COSPHI =DCOS( PHI  1 ) 
IF(DARS(COSPHI ).GT.1.0D0)  GO  TO  1 
RETURN 

1  continue 
call  flag(2, 6100, £200) 
100  continue 
20c  continue 

valuf  assigned  is  intfrim  value  only 
cosphi=1.0do 
kase  =  <; 

RETURN 
END 


FUNCTION  SINPHI(X) 

IMPLICIT  PEAL*8  (A-H,0-Z) 

REAL*B  K1,K2,K3,L 

SINPHI=DS0RT(1.0D0-(C0SPHI(X) )**2) 

PETUPN 

END 


FUNCTION    R(X) 
IMPLICIT    REAL*P     <A-H,0-Z) 
RFAI*P     K  1     K"?    K3 

C0MM0N/B1/E,H1,H2,P0IS,RAD1,RAD2,L,PHI 1 , PH I  2 , ALPHA , 
#THETA,KASE 
C0MM0N/R2/CPHI 1,CPHI2,C0,C1,C2 
IF(KASE.E0,2)   GO  TO  1 
R      =RAD14-L*(C0*X+(C1*X**2  )/2.000+(C2*X**3  )/3.0D0) 


IF(KASE.EQ.l) 
RFTURN 
1  CONTINUE 
VALUF  ASSIGNED  IS 
R=1.0D0 
RFTURN 
END 


R=RAD1 


INTERIM  VALUE  ONLY 


FUNCTION    G(X) 

IMPLICIT    REAL*8     {A-H,0-7) 

CCMM0N/B1/E,H1,H2,P0IS,RAD1,RAD?,L,PHU,PHI2,ALPHA, 
HFTA,KASF 


#THFTA,KASF 
Q=1.0D0/R(X) 
IF(KASE.EQ.2) 


Q=O.ODO 
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RETURN 
END 


FUNCTION    QC(X) 
IMPLICIT    PF4L*8     <A-H,0-Z) 
REAL*8    K! ,K2,K?,L 

COMMON /PI /P,H1 ,H2,POIStRADl, RAD?, Lt PHI  1 ,PHI2» ALPHA, 
#THFTA,KASF 
QO=GtX»/M  PHA 
IFIKASF.FQ.2)   OQ=THETA 
RFTURN 
FND 


FUNCTION  RKX) 

IMPLICIT  REAL*8  (A-H,0-Z) 

RFAL*8  K1,K2,K3,L 

COMMON/R!/E,Hl,H2,POIS,RADl,RAn2,L,PHIl , PHI  2, ALPHA, 
#THFTA,KASF 

COMMCN/B2/CPHI1 ,C PHI  2, CO, CI ,C2 

IF(KASF.GE.l)   GO  TQ  1 

R1=-L*SINPHI(X)/ (Cl+2.0DO*X*C2) 
C   THE  FOLLOWING  CARD  PROTECTS  AGAINST  OVERFLOW  IF  X  IS 
C   EXACTLY  AT  A  POINT  OF  INFLECTION 

IF ( PI .CT.1.0C50)  GO  TO  1 

RFTURN 
1  CONTINUF 
C   VALUF  ASSIGNED 

Ri=i.ono 

RETURN 
END 


IS  INTERIM  VALUE  ONLY 


FUNCTION    QKX) 
IMPLICIT    RFAL*8     <A-H,0-Z) 
REAL*8     Kl,K2,IO,L 

COMMON/B1/E, H1,H2,P0IS,RAD1,RAD2,L,PHI 1 , PHI  2, ALPHA, 
#THETA,KASE 
Q1=1,0DC/P1(X» 
IFIKASF.GF.l )   Ql=O.ODO 
RETURN 
END 

FUNCTION  DR(X) 

IMPLICIT  REAL*8  <A-H,0-Z) 

p  f.   A  t  *  P   |(  1   \SO      I/O   I 

CdMMnN/Pl/E,Hl,H2,POIS,RADl,RAD?,L,DHIl ,t>HT2,ALPHA, 
#THETA,KASE 
DR=L*C0SPHI( X) 
RETURN 
END 


1 


FUNCTION  DRl(X) 

IMPLICIT  REAL*8  (A-H, O-Z) 

COMMON/ 8 1 /E , Hit H2,P0IS,RA0t, RAD  2 ,L, PHI  1 , PHI  2, ALPHA, 
#THETA,KASE 
COMMON /B 2 /C PHI 1  ,CPHI2,C0,C1,C2 
IF(KASE.GE.l )   GO  TO  1 

DR1=DR(X) /SINPHI (X)+2,0D0*C2*(-R1(X ) )/ < C 1+2. 000* X *C 2 ) 
RETIJPN 
CONTINUE 
DR1=0.CDC 
RFTURN 
END 
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SUBROUTINE 
C   ELSTF  IS  THE 
C   ELEMENT  STIFP 
IMPLICIT  R 
REAL*8  Kl, 
DIMENSION 
#STIFEl(2At 
#STIFFF(24, 
COMMON/B1/ 
#THETA, KASE 
COMMCN/B2/ 
CPHIl=nCOS 
CPHI?=DCOS 
CO=CPHI1 

ci=2.onc*( 

C2=3.0DO*( 
IF(KASE.EQ 
C1=0.000 

C2=o.ono 

20  CONTINUE 
DO  50  1=1, 

no  so  j* if 

50  STK(I,J)=0 
C  ZEROS  OF  GAUSS 
XX1=. 96028 
XX2=.  79666 
XX3=. 52553 
XX4  =  . 18343 
C  COEFFICIENTS 
AA1=. 10122 
AA2=. 22238 
AA3=. 31370 
AA4  =  . 36268 
C  XG  CONTAINS  A 
C  COEFFICIENTS, 
WT(1 >=AA1* 
WT(2)=AA2* 
WT(3)=AA3* 
WT(4)=AA4* 
WT( 5)=WT<4 
WT(6)=WT<3 
WT(7)=WT(2 
WT(8)=WT(  1 
A=0.500 
XG( 1)=A+XX 
XG(  2)=A+XX 
XG(3)=A+XX 
XG(4)=A+XX 
XG(5)=A-XX 
XG(6)=A-XX 
XG(7)=A-XX 
XG(8)=A-XX 
DO  100  1=1 
00  100  J=l 
100  WTW(I,J)=W 
DO  200  1=1 
X=XG(  I  ) 
DO  200  J=l 
Y=XG( J) 
CALL  FORMK 
RALF=1.0C0 
DO  200  K=l 
DO  200  M=l 
200  STK(KtM)=S 
DO  300  1=1 
DO  300  J=l 
STKU  ,  J)  =  ( 
300  STK(J,I)=S 


XG( 8) ,WTW(3, 
,WT<  8), 


) 


ELSTF(STFMOD) 
CONTROL  PROGRAM  FOR  CONSTRUCTING  THE 
NFSS  MATRIX. 
EAL*8  (A-H,0-Z> 
K2  K3  i 

TINVT(24f24),TINV(24,24), 
?4)tSTF(24f  24),STK(24t  24) 
?4)  ,STFP0D(24,24) 
E,H1,H2,P0IS,RAD1,RAD2»L,PHI1 , PHI  2, ALPHA, 

CPHI1 ,CPHI2,C0,C1,C2 

(PHI1) 

(PHI2) 

3.CD0*(RA0  2-RAD1 >/L    -  2.  0D0*C PH 1 1 -C PHI  2  ) 
-2.0  00* (  PAD2-RAD1J/L  +  CPHIH-CPHT2  ) 
.01    GO    TO    20 


1*C.5D0 
2*C. 5D0 
3*0.500 
4*C.5D0 
4*0.500 
3*0.500 
2*C. 500 
1*0.500 

t8 

T( I)*WT( J) 

,8 


4 
4 

ODO 
AN  8 
8564 
4774 
4C99 
6424 
F  ZE 
5  362 
0344 
6458 
7833 
L  TH 
ALL 
.500 
.  500 
.  5D0 
.500 


POINT  FORMULAS  ARE 
9753623D0 
136267*00 
16328O9D0 
956498000 
ROS  ARE 
Q03762600 
533744700 
7788729D0 
783619800 
E  8  ZEROS  AND  WT  CONTAINS  ALL 

APE  SHIFTED  TH  RANGE  0  TO  1. 


THE 


,8 

(STF, X,Y) 

/QQ(X) 

t?4 

TK(K,M1+WTW<J,I)*STF(K,M>*L*RALF 

,24 

»24 

STK(  I,  JI+STM  J,  I  )  )*0.500 

TK(I ,J) 
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600 


IF(KASF.FQ.O)  RFTURN 

CALL  TFCRM(TINV,TINVT ) 

CALL    MUIT( TINVT,STK,STIFF1 |24,24t24) 

CALL    MULT(STIFF1,TINV,STIFFF,?A,24,2  4) 

00  600  T=l,24 

nn  60C  J=l,?4 

STIFFF(  I,  J)=0.5D0*(STIFFF(  I,J)+STIFFF(J,I)) 
STIFFF( J, I)=STIFFF( I, J) 

CONTINUF 

IFIKASE.EC?)     RFTURN 

CALL  RIGID  (STIFFF,STFMnD» 

RFTURN 

END 


SUBROUTINE  FORMK(STF,X,Y) 
FCRMK  FORMS  THE  ELEMENT  STIFFNESS  MATRIX. 
IMPLICIT  REALMS  (A-H,0-Z) 
REAL*8  K1,K2,K3,L 

DIMENSION  STF( 24,24)  , WORK1 ( 6 , 24  )  ,DFL ( 6 , 6 ) , 
CALL  DEL1(DEL,X) 
CALL  PST(P,X,Y) 


30C 


P( 6,24) 


DO  200 

on  200 

WORKK  I 
DO  200 


1=1,6 
J=l,?4 

,J)=O.ODO 
K  =  l,6 


200  WORKK I ,J)=WCPK1 (I,J)+DEL(I,K)*P(K,J) 


00 
DO 


^00 
300 


1=1, ?4 
J=l,24 


STF( I , J)=0.0D0 
00  ?00  K=l,6 
STF(  I, J)  =  STF(I 

RETUPN 
END 


J)+P(K, I )*W0RK1(K, J) 


350 


SURROU 
TFORM  FOR 
IMPLIC 
REAL*8 
DIMFMS 
COMMON 
#THETA, 
INITIALIZ 
DO  350 
DO  350 
TINVTf 
IF(KAS 
IF(KAS 
IF(KAS 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVTI 
TINVT( 
TINVT( 
TINVH 
TINVT( 
TINVT( 
TINVT( 
TINVTI 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVT( 
TINVT( 


TINE  TFORM(TINV,TINVT) 

MS  THE  TRANSFORMATION  MATRIX  AND  ITS  INVERSE, 

IT  REAL*8  (A-H,0-Z) 

K 1  K?  K"^  1 
ION fTINVT( 24,2  4) ,TINV<  24, ?4) 

/PI /E,H1,H2,P0IS,RAD1,RAD2,L,PHI1,PHI2, ALPHA, 
KASE 
E  TINVT  MATP  IX 

1=1,24 

J=l,?4 
I, J)=O.ODO 
E.EQ.2)  ALPHA=WIDTH 
E.EQ.2)  RAD1=1.0D0 
E.EQ.2)  RA02=1.000 


1 

1, 

I  , 

1» 

2, 

2, 

2, 

2, 

3, 
■x 

3 

n 

? 

■3 

4 

A 


1)=1. 
?)=-l 
3)=-l 
4)=1. 
5)  =  1 


ODO 
.  ODO 
.000 
ODO 
ODO 


6)=-1.0D0 
,  7)=-1.000 
,  8)=1.0D0 
,  9)=4.0D0 
,10)=-6.0D0 
,12)=2.000 
,13) =-6. ODO 
,14)=9.0D0 
♦16)=-3.0D0 
,21)=2.0D0 
,2  2)  =-3.  ODO 
,24)=1.0D0 
,  P)=2.0D0 
,  10)=-3.OD0 
4,12)=1.0D0 
4,13)=-A.0D0 
4,14)=6.0D0 


*L 
*L 
*L 
*L 
*L 
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TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TTNVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 
TINVT 


(   4, 

16) 

(      4, 

,17) 

(     4, 

IP) 

(     4  I 

,20) 

(      *>  i 

,    9) 

(      «=,  , 

,10) 

(      c;  1 

rill 

(      5, 

11) 

(    5, 

1*1 

(     5, 

15) 

(     5, 

21) 

(     5, 

??  ) 

(     c 

,231 

<     61 

9) 

(     6, 

,10) 

(     6, 

,11) 

(      A  , 

n) 

(    6, 

,14) 

(     6, 

,15) 

(     6, 

17) 

(     6, 

,1P) 

(     6, 

1^) 

(     7, 

I) 

(     7, 

,     ?) 

(     P, 

5) 

(     Pi 

,    6) 

(      Q, 

,     9) 

(     9, 

10) 

<    9, 

,12) 

(     9 1 

,13) 

(     ° 

,14) 

(     Q  > 

,16) 

(  10, 

,     9) 

(10, 

,10) 

(  10, 

,12) 

(  10, 

13) 

(  10, 

14) 

(10, 

,16) 

( 11 1 

9) 

(11 1 

,10) 

Ml, 

,11) 

(Hi 

13) 

<  11  1 

,14) 

(  1 1  1 

,15) 

(  12, 

,     9) 

(12, 

,10) 

(  I?, 

,11) 

(12, 

13) 

(12, 

,14) 

(  12, 

15) 

(1?, 

1  ) 

(  1*1 

,    5) 

(  15, 

9) 

(15, 

10) 

(  15, 

13) 

(  15, 

14) 

(16, 

9) 

(  16, 

,10) 

(  16, 

13) 

(  16, 

14) 

(  IT, 

,    9) 

(17, 

10) 

(  17, 

13) 

(  17, 

14) 

(IP, 

9) 

(IP, 

,10) 

(  IP, 

13) 

(  IP, 

,14) 

(  19, 

,     1) 

(1% 

3) 

(20, 

5) 

(20, 

,     7) 

=-2.0D0 

=  2.000 

=-3.000 

=  1.000 

=  2.000 

=-4.000 

=  2.000 

=  -3.000 

=  6.000 

=-3.0D0 

=1.000 

=-2.000 

=  1.000 

=  1.000 

=-2.000 

=  1.  000 

=-2.000 

=4.000 

=-2.000 

=  1.0  00 

=-2.000 

=  1.  000 

=  -1.000 

=  1.000 

=-1.000 

=  1.000 

=-4.000 

=6.000 

=-2.000 

=  6.0D0 

=-9.0D0 

=  3.000 

=  2.000 

=-3.000 

=  1.000 

=-2.000 

=  3.000 

=-1.000 

=-2.0D0 

=4.0  00 

=-2.000 

=3.000 

=-6.000 

=  3.000 

=  1.000 

=-2.000 

=  1.000 

=-1.000 

=2.000 

=-1.0D0 

=  1.000 

=  1.000 

=  4.0D0 

=  -6.000 

=-6.000 

=9.000 

=-2.0D0 

=  3.000 

=  2.000 

=-3.000 

=-2.000 

=  2.000 

=3.0D0 

=-3.0D0 

=  1.000 

=-1.000 

=-1.000 

=  1.000 

=-1.000 

=1.000 

=-1.000 

=  1.0D0 


*L 
*L 

*L 

*L 

*A!_PHA* 

*ALPHA* 

*ALPHA* 

*ALPHA* 

*ALPHA* 

*ALPHA* 

*ALPHA* 

*ALPHA* 

*ALPHA* 

*L*ALPH 

*L*ALPH 

*L*AIPH 

*L*ALPH 

*L*ALPH 

*L*ALPH 

*L  *ALPH 

*L*ALPH 

*L*ALPH 


RAD1 

RA01 

RA01 

RAD1 

RAD1 

RA01 

RAD! 

RA01 

RA01 

A*RAD1 

A*RA01 

A*RA01 

A*RA01 

A*RAC1 

A*RAO! 

A*RADl 

A*RAD1 

A*RAD! 


*L 

*L 

*L 

*L 

*L 

*L 

*AL 

*AL 

*AL 

*AL 

*AL 

*AL 

*l  * 

*\_+- 

*(_* 


PHA*R 
PHA*R 
PHA*R 
PHA*R 
PHA*R 
PHA*R 
ALPHA 
ALPHA 
ALPHA 
ALPHA 
ALPHA 
AL"H/\ 


AD2 

AD2 

AD2 

A02 

AD2 

AO? 

*PA02 

*RA02 

*RAD2 

*RA02 

*RAC2 

*RAD2 


*L 

*L 

*L 

*L 

*AL 

*AL 

*AL 

*AL 

*L* 

*L* 


PHA*RA02 

PHA*RA02 

PHA*RA02 

PHA*RAD2 

ALPHA*RAD2 

ALPHA*RA02 

ALPHA*RA02 

ALPHA*RA02 
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375 


SCO 


TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
TI 
DO 
00 
TI 
DO 
DO 
TI 
RF 
EN 


NVT( 
NVT( 
NVT( 
NVT( 
NVT( 
NVT( 
NVT  ( 
NVT( 
NVM 
NVT  ( 
NVT< 
NVT( 
NVT( 
NVT( 
NVT( 
NVT( 
MVT( 
NVT( 
NVT( 
NVT< 
NVT( 
NVTf 
NVT( 
NVT( 
375 
375 
NVT( 
500 
500 
NV(  I 

turn 

D 


21 
?! 
21 
21 
21 
2\ 
2  2 
22 
22 
?? 
22 
?2 
?3 
23 
2^ 
23 
23 
23 
74 
24 
24 
24 
24 
2** 


1  = 
J  = 
,  J 
1  = 
J  = 
J) 


q)=- 
10)=6 
13)=6 
14)=- 
21)=- 
2  2  )  =  3 

10)=3 
13)=4 
14)=- 
17)=- 
18)  =  3 

9)  =  2 
10)=- 
13)=- 
!  4 )  =3 
21)=1 
22)=- 

9)  =1 
10)=- 
13)=- 
14)  =? 
17)=1 
13)=- 
4,22, 
1,24 
)=-TI 
1,24 
It  24 
=  TINV 


4.000 

.000 

.000 

o.ODO 

2.0DO 

.000 

2.0D0 

.  ODO 

.000 

6.000 

2. ODO 

.000 

.000 

2.000 

3.000 

.000 

.000 

1.000 

.0  00 

1.000 

2.000 

,  ono 

.000 

1.000 

6 


*L 

*L 

*L 

*L 

*L 

*L 

*AL 

*AL 

*AL 

*AL 

*AL 

*AL 

*L* 

*  \_  A 

*|_  * 
■*L  * 


PHA* 
PHA* 
PHA* 
PHA* 
PHA* 
PHA* 
ALPH 
ALPH 
ALPH 
ALPH 
ALPH 
ALPH 


RAD1 

RAO! 

RAD1 

RAD1 

RADi 

RAD1 

A*RA01 

A*RAD! 

A*RAOl 

A*RAC1 

A*RAD' 

A*RAD1 


NVT(  I, J) 


T( J,  I  ) 


SUBROUTINE  DEL1(DEL,X) 
DELI  FORMS  THE  CONSTITUTIVE  LAW  MATRIX. 
IMPLICIT  REALMS  (A-M,n-Z) 
REAL*3  K1,K2,K3,L 
DIMENSION  DEL(6,6) 

C0MM0N/R1/E,H1,H2,P0IS,RAD1,RAD2,L,PHI  1 , PH I  2 , ALPHA , 
#THETAf KASE 
H( X)=Hl+( H2-H1 )*X 
Kl=(H(X)+F)/(!.0D0-POISA*2) 
K2=K1*PCIS 

K3=(K1*( l.ODC-POIS) ) /2.0D0 
01  =  (H(X)**3*E)/ (  1.2D1*( 1 .ODO-PO  IS**7  )  ) 

n->=m*pri  $ 

D3=(D1*< l.OHC-POIS) ) /2.0D0 
DO    ICO     1=1,6 
DO    100     J  =  l  ,6 
100    DEI ( I, J)=C.OOO 
0EL(1,1  )=K1 
OEL( 1 ,2)=K? 
DEL( 2,  1)  =  K7 
DEL(2,2)=K1 
DEL( 3,3)=K3 
DEL(4t4)=Dl 
DEL(4,5)=D2 
DEL( 5,4)=0? 
DEL(^,5)=Dl 
DEL(6,6)=D3 
RETURN 
END 

FUNCTION    P    IGKXtYl 
IMPLICIT    REAL**     (A-H,0-Z) 
R  ^  At  *■  R     Kl     K2     K3     I 

COMMON /R I Je, HI ,H2,P0IS,RAD1 ,RAD2,L,PHI1 , PH I  2 , ALPHA, 
4THETA,KASE 
P101  =Y/L 

RETURN 
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ENTRY 

P10? 

RETURN 

ENTRY 

P109 

RETURN 

FNTRY 

PllO 

RETURN 

ENTRY 

Pill 

RETURN 

ENTRY 

P112 

RETURN 

ENTRY 

P113 

RETURN 

ENTRY 

PI  14 

RETURN 

ENTRY 

P115 

RETURN 

ENTRY 

P116 

RETURN 

ENTRY 

P117 

RETURN 

ENTRY 

P118 

RETURN 

ENTRY 

Pll<? 

RETURN 

ENTRY 

P120 

RETURN 

ENTRY 

P121 

RETURN 

ENTRY 

P122 

RETURN 

ENTRY 

P123 

RETURN 

ENTdy 

P124 

RETURN 

ENO 


P    1G2(X,Y) 
=1.0DO/L 

P    109(X,Y> 

=Q1(X        )*X**3*Y**3 

P    110(X,Y) 

=Q1(X        )*X**3*Y**2 

P    U1(X,Y) 
=Q1(X        )*X**3*Y 

P    li2(X,Y| 
=G1(X       )*X**3 

P    U3(X,Y) 

=Q1 (X        )*x**2*Y**3 

P    llA(X.Y) 

=  G1(X        )*x**2*Y**2 

P    115(X.Y) 
=Q1(X        )*x**?*Y 

P     116(X,Y) 
=Q1(X        )*X**2 

P    117(X,Y) 
=G1(X        )*x*Y**3 

P    U8IX.YI 
=Q1(X        }*X*Y**2 

P    ll<MXtY) 
=Q1(X        )*X*Y 

P    120(X,Y) 
=G1(X       )*X 

P    121(X,Y) 

=  C? (X        )*Y**3 

P    122(X,Y) 
=G1(X       )*Y**? 

P    123(X,Y) 

=G1(X        >*Y 

P    12MX.Y) 
=G1(X        ) 


FUNCTION    P    2C1(X,Y> 
IMPLICIT    REAL*8     (A-H,0-Z) 
REAL*9    Kl ,K?,K3,L 
COMMON /Bl/E, HI ,H2 , PO I S , R ADl 
#THFTA,KASE 


RA02, L, PHI  1 i PHI?, ALPHA, 


P201 

RETURN 

ENTRY 

P202 

RETURN 

ENTRY 

P203 

RFTURN 

ENTRY 

P20A 

RETURN 

ENTRY 

P205 


=(CCSPHI(X        )*Q(X        ))*X*Y 


P    202(X,Y) 
=(CCSPHI(X 

P    203(X.Y) 
=(CCSPHI(X 

P    ?04<XfY) 
=(CCSPHI(X 


P    205( 
=  GO(X 


X.Y) 
)*X 


)  *  0  (  X        )  )  *  X 


)*Q(X        )I*Y 


)*Q(X        )) 
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RETURN 

ENTRY 

P207 

RETURN 

ENTRY 

P209 

RETURN 

ENTRY 

P210 

RETURN 

ENTRY 

P211 

RETURN 

ENTRY 

P212 

RETURN 

ENTRY 

P2!3 

RETURN 

ENTRY 

P2!4 

RETURN 

ENTRY 

°215 

RFTURN 

ENTRY 

P216 

RETURN 

ENTRY 

P?17 

RETURN 

ENTRY 

P21P 

RETURN 

ENTRY 

P21  °> 

RETURN 

ENTRY 

P?20 

RETURN 

ENTRY 

P221 

RETURN 

ENTRY 

P?2? 

RETURN 

ENTRY 

P223 

RETURN 

ENTRY 

P??4 

RETUPN 

END 


P    207( X,Y) 
=CQ(X        ) 

P    2C9(X,Y) 
=(SINPHI(X 

P    210(X,Y» 
=<SINPHI(X 

P    21KX,Y) 
=(SINPHI(X 

P    212( X,Y) 
=  (S  INPHKX 

P    213( X,Y) 
=  (SINPHKX 

P    214(X,Y) 
=  (S  INPHKX 

p     ?15(X,Y) 

=(SINPHI(X 

P    216(X,Y) 
=(SINPHI(X 

P     217(X,Y) 
=(SINPHI(X 

P    218( X,Y) 
=  (S  INPHKX 

P    ?19(X,Y) 
=(SINPHI(X 

P    220(X,Y) 
=  (SINPHKX 

P    22K  XtY) 
=  (SINPHKX 

P    222(X,Y) 
=(SINPHI(X 

P    223<X,Y) 
=  (SINPHKX 

P    224(X,Y) 
=  (SINPHKX 


)*Q(X 

)*0(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q<X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 
)*Q(X 


*  x**^*y**  3 
*X**3*Y*#2 
*X**3*Y 

*  X  -**  3 

*  X  *  *  2  *  Y  *  *  3 

*  X  *  *  2  *  Y  *  *  2 
*X**2VcY 
*X**2 
rX*Y**3 
-X*Y**2 
*X*Y 

*X 


FUNCTION    P    3  01 (XtY) 
IMPLICIT    REAL*8     (A-H,0-Z) 
R  F  Al  *  R     K  1     K?     K  3 

COMMON /Bl/E, HI , H2»POIStRADltRA02tL, PHI  1 ,PH 12 t ALPHA, 
*THETA,KASE 


P301 

=CQ(X   )*X 

RETURN 

ENTRY 

P  303( X,Y) 

P303 

=CQ(X   ) 

RETUPN 

ENTRY 

P  305(X,Y) 

°305 

Y/L-(Q(X 

)/L)*DR(X 

RETUPN 

ENTRY 

P  306C  X,YI 

P306 

=1.000/L-(0(X 

)/L)*OR(X 

RETURN 

ENTRY 

P  307 ( X,Y) 

)*X*Y 


)*X 
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P307 

RETURN 

ENTRY 

P308 

RFTURN 

END 


=(-C(X 


P    30R(XTY  ) 

=  (-C(x      »/d*or(x 


)/L)*DR(X        »*Y 
) 


FUNCTION    P    4C1 ( X,Y) 
IMPLICIT    REAL*8    (A-HtO-Z) 
R F  Al  ♦  8     K  1      K  ?     K  "^     I 

COMMON /Rl/E, HI  , H? t PO I S , R AD  1 » RAO? , L t °HI 1 , PH I  2 , ALPHA, 
*THETA,KASE 


P401 

RFTURN 

ENTRY 

P402 

RETURN 

ENTRY 

P403 

RETUPN 

ENTRY 

P404 

RFTUPN 

ENTRY 

P409 

RETUPN 

ENTRY 

P410 

RFTURN 

ENTRY 

P411 

RETURN 

FNTRY 

P412 

RETUPN 

FNTRY 

P413 

RETUPN 

ENTRY 

P414 

RETURN 

ENTRY 

P415 

RETURN 

ENTRY 

P416 

RFTUPN 

ENTRY 

P490 

RFTURN 

FNTRY 

PAQ1 

RFTUPN 

ENO 


=  -(Gl(X   )**2/L)*IWl(X 

P  40?( X»Y) 

=  -(Cl(X   )**2/L)*DPKX 

P  4C3<X,Y) 

=-(Cl(X   )+*2/L)*ORl(X 

o  404(X,Y) 

=-(Ql(X   )**2/L)*DRl(X 

P  40Q( X,Y) 
=<-6.0DO/L**2  )*X*Y**3 

p  410( X,Y) 
=(-6.0DO/L**2  )*X*Y**2 

P  411(X,Y) 

=(-6.000/L**2  )*X*Y 

P  412( X,Y) 

=  (-6.000/L**2  )*X 

P  413(X,Y) 
=(-2.OD0/L**2  )*Y**3 

P  414( XTY  ) 
=(-2.0DO/L**?  )*Y**2 

P  415(X,Y  ) 
=(-2.0DO/L**2  )*Y 

P  4  1 6  (  X  ,  Y  ) 
=(-2.OD0/L**2 ) 


P  490(X,Y) 
=COSPHI (X 

P  4<U(  XrY» 
=CQ(X   )**? 


)*0(X   )/L 


)*X *Y+(01( X 
)*X  *-(QUX 
)*Y 

) 


)/L»*Y 


FUNCTION  P  5C1(X,Y) 
IMPLICIT  REAL*R  <A-H,0-Z) 
RFAL*8  Kl  K2  K3  L 

COMMON/B1/E,H1 ,H2 , PO I S , R AOl , RAO? ,L , PHI  1 , PH T 2 , ALPHA, 
#THETA,KASE 

)*Q(X   )*Q1(X 


P501 

RETURN 

ENTRY 

P5D? 

RETURN 

ENTRY 

"503 

RETUPN 

ENTRY 

P504 


=COSPHI(X 

P  S02(XtY» 
=COSPHI (X 

P  503(X,Y) 
=  C0SPHI  (X 

P  504(X,Y) 
=COSPHI (X 


)+Q(X 
)*Q(X 
)*Q(X 


)*Q1(X 
)*Q1 (X 
)*01(X 


)*X*Y 
)  "X 
)*Y 
) 
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RETURN 
FNTRY 
P505 
RETURN 
ENTRY 
P507 
RETURN 
ENTRY 
P509 
*X**-**Y 
RETURN 
ENTRY 
P510 

RETURN 

ENTRY 

0511 

RETURN 

ENTRY 

P512 

RETURN 

ENTRY 

P513 

RETURN 

FNTRY 

P51A 

RETUPN 

ENTRY 

P515 

RETUPN 

FNTRY 

P516 

RETUPN 

ENTRY 

P517 

RETURN 

ENTPY 

P518 

RETURN 

ENTRY 

P519 

RETURN 

ENTRY 

P520 

RETURN 

ENTRY 

P521 

RETUPN 

ENTRY 

P522 

RETURN 

FNTRY 

P590 

RETURN 

ENTRY 

P5<H 

RETUPN 

ENTRY 

P592 

RETURN 

END 


P  505(XtY) 

=  Q(X   )*QQ(X   )*X 

P  507(X,Y) 
=G(X   )*GQ(X   ) 

P  509<X,Y) 

=  -3.0DO*P490(Xt  Y)*X**2*Y**3-6.  0D0*P491  (  X,  Y)* 


P  510(X,Y) 

=-3.0D0*P490( X, Y)*X**?*Y**2-2.0D0*P491  (X,Y  )* 


P  511 ( X,Y) 

=  -3.0nO*P49  0(X,Y)*X**'"tY 

P  512( X,YI 
=-3.0D0*P490(X, Y)*X**2 

P  513( X, Y) 

=-2.0D0*P49O<X,Y)  -*X*Y*  -«3-fe.  00  0*P4  91  (  X  ,  Y  )  *  X*  *  2 

P  514(  X,  Y  ) 

=-?.0D0*P49  0( X, Y) *X*Y**?-2.0D0*P491 (X,Y  J*X**2 

P  515(X,Y  ) 
=-2.0D0*P49  0( X,Y) *X*Y 

P  516( X,Y ) 
=-?.ODO*P490( X,Y) *X 

P  517(X,Y) 

=  -P4O0(X,Y)-*Y**3-5.0Q0*P4Ql(X,Y)*X*Y 

P  518(X,Y ) 

=-PA90( X,Y)*Y**2-2.0D0*P4°1(X»Y)*X 

P  519(X,Y) 
=-P490(X, Y) J  Y 

P  ?20(X,Y) 
=  -t>A90(X,  Y) 

P  521 (X,Y) 
=-6.0D0*P491  ( X,  Y)*Y 

P  ^22 (X,Y) 
=-2.0D0*P491 (X,Y) 

P  590(X,Y) 

=(Q(X   )**2/L)*HR(X   ) 

P  591 (X,Y) 

=(0(X   )*QQ(X   )/L)*DR(X   ) 

P  592(X,Y) 
=QQ(X   )/L 


FUNCTION  P  6CK  X,Y) 
IMPLICIT  REAL*8  (A-HfO-Z) 
REAL*^  K1«K?  K3  L 

C0MM0N/B1/E,H1,H2,P0IS,RAD1,RAD2»L,PHI1 , PH I  2 » ALPHA, 
#THETAtKASE 
P601      =(QC(X   )*Q1<X   ))*X 
RETURN 
ENTRY     p  603<X,Y  ) 
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P603 

RETURN 

ENTRY 

P60^ 

RFTUPN 

EN^PY 

P606 

RFTUPN 

ENTRY 

P607 

RETURN 

CNTPY 

P608 

RETURN 

ENTRY 

#X**2*Y 
RETURN- 
ENTRY 
P610 

#X**2*Y 
RETURN 
ENTRY 
P611 
RETURN 
ENTRY 
P613 

«X*Y**2 
RETURN 
ENTRY 
P614 
RETURN 
ENTRY 
P615 
RETURN 
FNTRY 
P617 
RETURN 
ENTRY 
P61R 
RETURN 
ENTRY 
P619 
RFTurn 

ENTRY 

P621 

RETURN 

ENTRY 

P622 

RETURN 

ENTRY 

P623 

RETURN 

END 


**2 


=  CO(X       1  +  QHX       ) 

P    605( X,YI 

=  -2.0D0*P590(X,Y)*X*YMQ(X        )/L)*Y 

P    6  06 (X,Y) 

=  -?,0r>0*P5QrMXtY)  *X       +(0(X       )/L) 

P    607(X,Y) 

=    -2.000*P590(X,Y)*Y 

P    608(X,Y) 

=    -2.000*P590(X,Y) 

P    609(X,Y) 

=  6. 0DO*P5<n  <XtY)*X**^*Y**2-18.0D0*P59? (X, Y)* 


P    610( X,Y) 

=4.0D0*P5  91 (X,Y)*X**3*Y         -12- 0P0*P592( X, Y)* 


P    61i(  X,Y) 

=2. 000* P5 91 (X, Y)*X**3-6.0D0*P5  92<  X, Y»*X**2 

P    613( X,Y) 

=6.0D0*P5  91  (X,Y  )*X**2*Y**2-12.OD0*P592(  X,  Y)* 


P    614(X,Y) 

=  4.  CnO*P5  91  (X,Y)*X**2*Y-8.0D0*P59  2(X,Y)*X*Y 


P    615( X,Y) 

=2.0DO*P5  91(XtY)*X**2 


-    4. 0P0*P59? <X,Y) *X 


P    617(  X»Y) 

=6. CDO^PSQl ( X, Y)*X*Y*-  2-6. ODO*P5P2( X,Y) *Y**2 

P    618< X,Y) 

=4. 0D0*P591( X,Y)*X         *Y         -    4.  0D0*P59?  (  X  ,  Y  )  *Y 


P    619( X,Y) 
=2.CDO*P591( X,Y)*X 

P    621 ( X,Y) 

=6. COO* P 5 91 (X, Y)*Y**2 

P    622( X,YI 
=4.0D0*P591(X,Y)*Y 

P    6?3( X,Y) 
=2.0DO*P591(X,Y) 


-    2.0DO*P5Q2 (X»Y) 


SUBROUTINE  PST(P,X,YI 
C   PST  CONVERTS  THE  P»S  FROM  FUNCTIONS 
C   FOR  ELEMFN'TS  OF  THE  P  MATRIX, 

IMPLICIT    REAL*8     <A-HfO-Z) 

REAL*8  K1,K2,K3,L 

DIMENSION  P(6,24) 


OF  X  AND  Y  TO  VALUES 


P(l 
P(l 
P(  1 
P(l 
P<1 
P(  1 
P(l 
P(l 
P(l 
P(l 


,  n=pioi(x,Y) 

,  2)=P102(X,Y) 

,  3)=0.0DC 

,  41=0. ODO 

,  5I=C.CDC 

,  6)=C.0DC 

,  7)=C.0D0 

,  P)=O.0DC 

,  9)=P1G9(X,Y) 
»10»=P110(X,Y) 
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RETURN 
FNTRY 
P505 
RETURN 
ENTRY 
P507 
RETURN 
ENTRY 
P509 
#X**"**Y 
RETURN 
ENTRY 
P510 

RETURN 
PNTRY 
P511 
RETURN 
FNTRY 
P512 
RFTURN 
ENTRY 
P513 
#*Y 
RETURN 
FNTRY 
P514 
RFTURN 
ENTRY 
P515 
RETURN 
FNTRY 
P516 
RETURN 
ENTRY 
P517 
RETURN 
ENTRY 
P518 
RETURN 
ENTRY 
P519 
RETURN 
ENTRY 
P520 
RFTURN 
ENTRY 
P521 
RETURN 
ENTRY 
P522 
RETURN 
FNTRY 
P59Q 
RETURN 
FNTRY 
P591 
RETURN 
ENTRY 
P592 
RETURN 
END 


P  505( X,Y) 

=G(X   )*QQ(X   )*X 

P  507<X,Y) 
=G(X   )*00(X   ) 

P  509<X,Y) 

=  -3.0D0*P490(X,Y)*X**2*Y**3-<S.0D0*P491(  X,  Y)* 


P  510(X,Y) 

=-^.00  0*P490( X, Y)*X**2*Y**2-2.000*P491  (X»Y  )* 


P  51K  X,Y) 

=  -3.0nO*P49  0(X,Y)*X**">*Y 

P  512< X,Y) 

=-3.0D0*P490(X,Y)*X**2 

P  5131 X,Y) 

=-2.0nO*P490(X,Y) *X*Y**?-6. 000*P491(X,YI*X**2 


P  514( X,Y  ) 

=-?.ODO*P490(X, Y) *X*Y**2-2.0D0*P491 (X, Y )*X**2 

P  515(X,Y  ) 
=-2.0D0*P490< X, Y) *X*Y 

P  516( X,Y  ) 
=-?.0D0*P490( X,Y) *X 

P  517(X,Y) 
=-PAO0(X,Y)«Y**3-^.0D0^P4Ql(X,Y)*X*Y 

P  518(X,Y  ) 

=-P490(X, Y)*Y**?-2.0D0*P4Q1 (X,Y)*X 

P  519(XfY) 
=-P490(X, Y) J   Y 

P  520( XrY  ) 
=  -c>490(  X,  Y) 

P  521 (X,Y) 
=-6.0D0*P491 ( X, Y) *Y 

P  522(X,Y) 
=-2.000*P491 (X, Y) 

P  590(X,Y) 

=(Q(X   )**2/L)*HP(X   ) 

P  591 (X,Y) 

=(0(X   )*QQ(X   )/L)*nR(X   ) 

P  592(X,Y) 

=QQ(X   )/L 


FUNCTION  P  6C1( X,Y) 
IMPLICIT  REAL**  (A-H,0-Z) 
REAL**  K1  K?  K^  I 

C0MM0N/B1/E,H1,H2,P0IS,RAD1,RAD2,L,PHI1 , PHI2 , ALPHA, 
#THETA,KASE 
P601      =(QQ(X   )*Q1(X   ))*X 
RETURN 
ENTRY     p  603<X,Y) 
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P603 

RETURN 

ENTRY 

P60* 

RFTUPN 

FNTRY 

P606 

RFTUPN 

ENTRY 

P607 

RETURN 

CNTPY 

P608 

RETURN 

FNTRY 

P609 

#x+*2*y**; 

RETURN 

ENTRY 

P610 
#X**2*Y 

RETURN 

FNTRY 

P611 

RETURN 

ENTRY 

P613 
*X*Y**2 

RFTUPN 

ENTRY 

P614 

RFTURN 

ENTRY 

P615 

RFTUPN 

FNTRY 

P617 

RETURN 

ENTRY 

P618 

RETURN 

ENTRY 

P619 

RFTURN 

ENTRY 

P621 

RETURN 

ENTRY 

P622 

RETURN 

ENTRY 

P623 

RFTUPN 

END 


=  CO(X   ItQKX   ) 

P  605( X,Y» 

=-2.0D0*P59  0( X,Y) *X*Y  + (Q(X   ) /L)*Y 

P  6  06( X,Y) 

=  -?.ODO*P5QO(X,Y)  *X   +(Q(X   )/L) 

P  607( X»Y) 

=  -2.0D0*P590(X,Y)*Y 

P  608< X,Y  > 

=  -2.0DO*P590(X,Y> 

P  609(X,Y) 

=  6. 0D0*P5<U  (XtY)*X**^*Y**2-18.0DO*P592 (X, Y)* 


P  610(  X,Y) 

=4.0D0*P5  91(X,Y )*X**3*Y    -  12. ODO*P 592 ( X,Y)* 


P  6  1 1  <  X  ,  Y  » 

=2.0D0*P5  91 <X»Y)*X**3-6.000*P592< X, Y)*X**2 

P  613( X,Y) 

=6.0D0*P591  (X,Y  »*X«*2*Y**2-12.0DO*P592(  X,  Y»* 


P  614( X,Y) 

=  4. CD0*P5  91 (X,Y)*X**2*Y-8.0DO*P59?(X,Y )*X*Y 


P  615( X,Y) 

=  2.  0DO*P5  91 ( X,Y)*X**2 


-  4. 0P0*P59? (X,Y)*X 


p  617<  X»Y) 

=  6.  CDO*P5<n  (  X,  Y)*X*Y*'  2-6.  ODO*P5^2(  XfY)*Y**-2 

P  618(  X,Y) 

=4. 0D0*P591( X,Y)*X    ^Y    -  4.  ODD*  P59?  (  X  ,  Y  )  *Y 


P  619( X,Y) 
=2.CDO*P591( X,Y)*X 

P  621 ( XtY) 

=6. COO* P 591 (X, Y)*Y**2 

P  622( X,Y> 
=4.0D0*P591(X,Y)*Y 

P  6?3( X,Y) 

=2.0DO*P591(X,Y) 


-  2.0DO*P592<XfY) 


SUBROUTINE  PST(P,X,Y) 
C   PST  CONVFRTS  THE  P«S  FROM  FUNCTIONS  OF  X  AND  Y  TO  VALUES 
C   FOR  ELEMENTS  OF  THE  P  MATRIX. 

IMPLICIT    REAL*8     (A-H,0-Z) 

RFAL*8  KlfK2,K3,L 

DIMENSION  P<6,24) 

p(i,  n  =  Pini(x,Y) 

P(lt  2>=P102(X,Y) 

PCI,  3)=0.0DC 

P(l»  41=0. ODO 

P(lt  5»=C.CDC 

P(l,  6)=C.0DC 

Plli  71=0. ODO 

P(lt  8)=O.0DC 

P(l,  9)=P1G9(X,Y) 
P(l  ,10»  =  P110(X,Y) 
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P 
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P 
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P 

O 

c 
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P( 

2, 

22) 

P 

2, 
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PI 

2, 
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PI 

3, 

1) 

PI 

3, 

2) 

D  1 

3, 

,     3  ) 

PI 

3, 

A  ) 

PI 
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,     ^  ) 

PI 

3, 

,     6) 

PI 

3, 

7) 

PI 

^  , 

8) 

PI 

3, 

,     9) 

PI 

3, 

10) 

PI 

3, 

11  ) 

PI 

3, 

12) 

PI 

3, 

13) 

PI 

3  , 

14) 

P( 

3, 

15) 

PI 

3, 

16) 

P( 

3, 

17) 

PI 

3, 

18) 

PI 

3, 

,1<?) 

PI 

3, 

20) 

PI 

3, 

21  ) 

PI 

?, 

,22) 

PI 

3, 

23) 

PI 

3, 

24) 

PI 

4, 

1  ) 

P( 

A  , 

2) 

PI 

4i 

3  ) 

PI 

4, 

4) 

P( 

4, 

5) 

PI 

^  i 

6) 

PI 

4, 

7) 

P( 

A, 

8) 

P< 

A, 

<3) 

PI 

4, 

10) 

=P1 11 (X,Y) 
=P! 12(X,Y) 

=pinix, 

=P114(X, 
=  P11MX, 
=  P  1  1  6  ( X  , 
=P!17(X, 

=pinix, 

=P119(X, 
=P1?0(X, 
=  P 1  2 1  (  X  , 
=  P122 (X, 
=P123(X, 
=P124(X, 
=P201 (X, 
=P2C2(X, 
=  P20^(X, Y  ) 
=  P?CMX,  Y> 
=P205(X, Y) 

=c.ono 

=P2C7(X, Y) 

=0. one 

=P203(X, 

=P210(X, 

=P211(X, 

=P212 (X, 

=P213(X, 

=P214(X, 

=P215(XT 

=P216(X, 

=P217(X, 

=P218(X, 

=  P21(3(X, 

=P220(X, 

=  P221 (X, 

=P222(X, 

=P223(X, 

=P??4(X, 

=P301(X, 

=  C000 

=P3C3(XrY) 

=C.0O0 

=P305(X 

=P^C6(Xt 

=P3C7( X, 

=P303(X, 

=0. ODO 

=O.CDO 

=0.0D0 

=C.CDC 

=0.000 

=0,000 

=  c.coo 

=0.000 
=O.0DO 

=c.coc 

=0.000 
=0.000 

=c.coc 

=0.000 
=C.CDC 
=C. CDC 
=  P401  (X 
=P402(X, Y) 
=PA03(X, Y) 
=P404(X, Y) 

=o.ono 
=o.ooo 

=0.000 

=c.ono 

=P409(X, Y) 
=P410(X,Y) 


Y) 
Y) 
Y) 
Y) 
Y) 
Y) 
Y) 
Y) 
Y) 
Y) 
Y) 
Y  ) 
Y) 
Y) 


Y  ) 

Y) 

Y) 

Y) 

Y) 

Y) 

Y 

Y 

Y 

Y 

Y 

Y 

Y) 

Y) 

Y) 

Y) 

Y) 


Y) 
Y) 
Y) 
Y) 


Y) 
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p 
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p 

|4 
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►     1) 

p 

S, 
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p 
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PI 
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p 
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i     5  ) 

p 
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►     7) 
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►  1C) 

p 

5, 

►  11  ) 
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5, 
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PI 

5, 

►  "i  "*  ) 

P. 

[5, 

►  14) 

PI 

5, 

►  15) 

P 

5, 

►  16) 

Pi 

5, 

►  17) 

PI 

S, 

►  18) 

PI 

5, 

,19  ) 

Pi 

5  - 

►  20) 

P< 

S, 

►  21) 

PI 

5 

►  22) 

p  i 

5, 

►  23) 

PI 

5, 

►  24) 

PI 

61 

►    1  ) 

P 

fti 

►    2) 

PI 

6, 

►     7  ) 

PI 

6, 

►     4  ) 

PI 

6, 

►     5) 

P( 

6, 

►     6) 

o| 

6, 

►    7) 

PI 

6, 

►     8) 

P( 

fti 

9) 

PI 

6, 

.10) 

PI 

6, 

.11) 

P( 

6, 

12) 

PI 

ft  1 

13) 

P( 

6, 

14) 

p  1 

fti 

15) 

PI 

6, 

1ft) 

P( 

6, 

17) 

P( 

fti 

18) 

P( 

ft  1 

.19) 

p( 

fti 

20) 

P! 

61 

21  ) 

P( 

fti 

22) 

P( 

61 

23) 

P( 

fti 

24) 

RETURN 
END 


=P411(X,Y) 

=P412(X,Y) 

=P413(X,Y) 

=P414(X,Y> 

=P415(X,Y) 

=P41ft(X,Y) 

=C.0D0 

=0,000 

=0.000 

=C.0DC 

=0.000 

=c.coo 
=c.ooc 

=C.0D0 

=P5C1(X,Y) 

=P502(X, Y) 

=P503(X,Y) 

=P504(X,Y) 

=P505(X»Y) 

=0.000 

=P5C7(X, Y) 

=0.000 

=P509(X,Y) 

=P510(X, Y) 

=P511(X,Y1 

=P512(X,Y) 

=P513(X, Y) 

=P514(X, Y) 

=P515(X,Y) 

=P51ft(X, Y) 

=P517<X, Y) 

=P518(X,Y) 

=P519(X, Y) 

=P520(X,Y) 

=P5?1(X,Y) 

=P522(X,Y) 

=0.000 

=0.000 

=P601 <X,Y) 

=0.000 

=P6C3(X, Y) 

=0.000 

=P6C5(X,Y) 

=PftC6(X, Y) 

=Pft07(X,Y) 

=Pft08(X,Y) 

=P609(  X,  Y) 

=PftlO(X, Y) 

=P611 <X,Y) 

=0.000 

=Pftl3(X,Y) 

=P614(X, Y) 

=Pftl5(X, Y) 

=0.000 

=P617( X, Y) 

=Pftl8(X,Y) 

=Pftl9(X, Y) 

=C.0O0 

=P621(X,Y) 

=Pft22 (X,Y) 

=P623(XtY) 

=0.000 


SUBROUTINE  RIGID  { ST  I FFF , STFMOD  ) 
RIGID  FORMS  THE  RIGID  BODY  AND  THE  K22(TKT)  MATRIX. 
IMPLICIT  REAL*8  (A-H.0-7.) 
R  F  A  t  ♦  8  Kl  K?  K^ 

DIMENSION* TTf( 24,6 ) ,STFRB( 24,24) »AKT (24 ,ft) tTKTINV(6,ft) 
#, RBI  24,6) , ST IFFF( 2  4,24 )f ST F MODI 24,24),TKT(6,6),AK(3ft), 
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#AKTT(6, 
COMMON/ 
«THFTA,K 
COMMON/ 
EQUIVAL 
F  Z  (  X  )  =  L 
CO  =  CPHI 
IMKASF 
IF(KAS^ 
IF(KASF 
IF(KASF 
IF(KASF 
SPHI1=D 
SPHI2=D 
NIJMPRICALL 
A=0. 50C 
C=.4R01 
7=.506l 
C=.^9*3 
7=7*.  11 
C  =  . 26?7 
Z  =  7  +  . ] 5 
C=.9171 
Z  =  Z  +  .  1Q 
TFfKASE 
IF(KASF 
IF(KASF 
INITIALIZF 
1  = 
J  = 
t  J 
1  = 
J  = 
J) 
1 


24) 

PI/ 
ASF 
B2/ 

FNC 

*ds 
l 

•  GE 

.Gr 

•  GF 

.FQ 

•  GE 

SQR 
SOR 
Y    F 


c,Hl,H2,eniS,RADl,RAD?,L, PHI  1  , PHI 2, ALPHA, 

CPHI1  ,CPHI2  ,CO,Cl,C2 
E     <TKT( 1,1 1  ,AK(1 ) ) 

CRT(  l.OOO-fCO  +  Cl^X-K^X**^)*^) 


.1)   ri=o.ooo 
.1 )  r.2  =  o.ono 

•  1  )    CPHIl  =  O.ODO 

.3)    CPHI 1=(RAD2-RA01 ) /L 

.1  >    CPHI2  =  CPHI1 

T(  L.0P0-CPHU**2) 

T(  I.OOO-CPHI?**?) 

VALUATE    FUNCTION    Z. 


no  i 
no   l 

TKT(  I 


00 

no 

RB 
RB 
RR 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RR 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 


2 

-> 

If 

1 


1 

1 

2 

2 

2 
3 

3 
3 
4 
c 

4 

5 
f 
6 

7, 

7, 

7, 

B, 

B, 

Q, 

°t 

°t 

10, 

10, 

10, 

11 

11 

12 

12 

13 

13 

13, 

13, 

13, 

14, 

14, 


3 

2 
4 
6 
1 
3 
5 
1 

? 

1 
^ 

5 

2 

6 
1 
3 

1 
3 
5 
t2 
,4 
»2 
t 


44 

42 

i  -j 

if 
6f 
^fi 

73 
13 
.E 
.E 

.F 

R 

It 
1  , 
)  = 
It 
It 
=  C 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 
)  = 


7  = 
Z  = 

BO 


q?q?4R7 
6R1 A51 P 
23R7068 
90^172? 

?  a v -581 

R  ? 3 22°  3 
2  1  ?  4  7  B  ? 
41«sl6fl 
0.1  )  Z  = 
0.2  > 
0.3) 
IGIC 
6 
6 

O.ODO 
?4 
6 

.000 
CPHI1 
-SPHI 1 
Z*RB( 
1.000 
-Z 
QAD1 
SPHI1 
RB<  1, 
Z  +  RB( 
-(  (  RB 

Cl/L 

-RB(3, 

RB(  3, 

-RB(3, 

RB(4, 

-RB( 4, 

CPHI2 

-SPHI2 

-RB(7t 

1  .000 

RAD2 

SPHI2 

CPHI2 

-RAD2* 

-( (RB( 

(CU2 

-RBP, 

RB(q, 

RB(9, 
RB(  10 

-RB(10 
R8(7, 
RB(7, 
RB(7, 
RAD2* 
"A02* 

-DSIN( 
OCOS( 


6P12no 

Rl  ROn-1  MFZ(  A+f  J+FZC  Kl) 

! 337no 

66R724n0*(FZ(A  +  r)+FZ(A-O) 

6449nO 

PQ43ft4nO*(FZ( A+C)+F7( A-C ) ) 

49on-i 

^lROc>c5nO*(F7(A  +  C)  +  FZ(A-C)) 

L 

L 

L*SPHI1 

OY  MATRIX  AND  K??(TKTI  MATRIX, 


1,  1  >-RADl*RB( 1, 3) 


1) 

3,  1  )-PA01*PB(3,  1) 

(1  ,1)/RB(3,1)  )*C1  )/L 

1  )*v2+Z*R8(4,  I  )-CO*RB(  1  ,  1  )  -  (  R  AD1  *C  1  )  /L 

I)  /RAD1 

5)/RA01 

l)/OAni 

^) /RA01 


3)*RA02 


RB 
7, 
.0 
1) 
1) 
^) 
,1 
,5 
1) 
1  ) 
3) 
OS 
RB 
AL 
AL 


(9,3) 

1  )/RB( 9,  1)  )*(CH-?.0D0*C2)  )  /L 

D0*C2)/L 

**2-<C0+Cl*C2)*RB(9,3)-(RAD2J'(CH-2.0 

/RAD2 

I/RAD2 
)/PA02 

*OCOS( ALPHA) 
*DSIN( ALPHA) 

IN(ALPHA)*RB(7»3I 

(9,1 )*DCOS( ALPHA) 

PHA) 

PHA) 
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RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RR 
RR 
RR 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RR 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
DO 

no 

RB 
IF 
IF 

DO 
DO 

JL=?+J 
JR=3+J 
RB(  I, J 
IF(KAS 
GO    TO 
NDIM=6 
CONTIN 
DO    10 
DO    10 
AKTU  , 
DO     10 
AKT(  I, 
CONTIN 
DO    20 
DO    20 
TKTU  , 
DO    20 
RB(K,I)=R 
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1* 
15 
15 
1^ 
15 
15 
16 
16 
16 
16 
16 
17 
17 
17 
17 
18 
1^ 
IB 
1R 
19 
19 
10 

1° 
19 
20 
20 
20 
20 
20 
21 
21 
21 
21 
21 
22 
22 
22 
22 
22 
23 
23 
23 
23 
24 
24 
24 
24 
3  I 

3  J 
I, J 
KAS 
KAS 

4  I 
4  J 


6)  = 
1  )  = 

2)  = 

3)  = 

4)  = 

5)  = 

1)  = 

2)  = 

3)  = 

4)  = 

5)  = 
1  )  = 
2)  = 
4)  = 
51  = 
1  )  = 
2)  = 
4)  = 
5>  = 

1)  = 

2)  = 

3)  = 
A)  = 
51  = 

1)  = 

2)  = 

4)  = 
c)  = 

)  = 
)  = 
)  = 
)  = 

4)  = 

5)  = 

1)  = 

2)  = 
?)  = 

4)  = 

5)  = 
1  )  = 
2)  = 

4)  = 

5)  = 
II* 
2)  = 

4)  = 

5)  = 
=  4, 
=  1, 
)=- 
E.E 
F.G 
=  lt 
=  lt 


RAD2 


RB 

-RR 
RB 
RB 

-RB 
RB 

-RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 
RB 

-RB 
RB 

-RB 
RB 
RB 
RB 
RB 
RB 
RB 

-RB 
RB 
RB 

-RB 
RB 
RB 
RB 
RB 

-RB 
RB 
RR 
RR 
RB 
RB 
RB 

-PR 

RB 

RB 

RB 

-RB 

22,6 

6 

RB(  If 

0.0) 

T.l  ) 

24 

3 


9, 
7, 
13 
13 
10 
10 
10 
14 
14 
9, 
15 
13 
13 
10 
16 
14 
14 

It 

3, 
It 
14 
14 
1  A 
14 
2, 
2, 
2, 
3, 
3, 
It 
14 
14 
4, 
4, 
4, 
14 
14 
5, 
5, 
14 
14 
4, 
22 
14 
14 


J) 

GO 

ND 


1)*RR( 

1)*RB( 

1) 

,  2  )  *  P  A 

,1)*RA 

,1  )*RB 

,1  )*RB 

,3) 

,1  )*RB 

,2)*RB 

1)*RB( 

,1  )/RA 

tl) 

t2) 

,1  )*RB 

,l)/RA 

,2)*RB 

,1  )*RB 

1)*RB( 

3)+RB< 

3) 

tl  )*PB 

,2)*Rfi 

tl) 

,2) 

4)*RB( 

4)*RB( 

6) 

1)*RB( 

i)*RB( 

1) 

♦  1  )*PB 

,2)*RB 

1  )*RB( 

1)*RB( 

3) 

,1  )*RB 

,2)*RB 

2) *RR( 

2 )*RB( 

,2)*RB 

t 1  )*PR 

1)*RB( 

,1)/RA 

,2)*PB 

,1)*RB 


14,2) 
14,1) 

D2 
D2 

(14,2) 
(14,1 ) 

(10,5) 
(10,5) 
14,1 )/RAD2 
D2 


(14,  D/RAD2 

D2 

(12,4) 

(12,4) 

14,2) 

14,1  ) 

(  1,5) 
(1,5) 


14,2) 
14,1) 

14,2) 
14,1) 

(3,5) 
(  3,5) 
14,2) 
14,1  ) 

(4,5) 

(4,5) 

14,1  ) 

14,2) 

(5,4) 

(5,4) 

14,1  )/RADl 

Dl 

(6,4) 

(6,4) 


TO  5 
IM=5 


L)  =  RB(  I, JR  ) 

E.EQ.l  )  NDIM=4 

7 


UE 
1  = 
J  = 
J) 
K  = 
J) 
UF 
1  = 
J  = 
J) 
K  = 
B( 


1,24 

1,NDIH 

=0.0D0 

lt24 

=  AKT(I, J)+STIFFF(  I,K)*RR(K,J ) 

1,NDIM 

1 ,NDIM 

=C.0D0 

lt?4 

I,K)  TRANSPOSE 
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20 


21 


22 

25 


50 


6C 

1000 


TKT 
COM 
IF( 

on 

1ST 

nn 
jp= 

AK( 
IF( 

on 

1ST 

no 
jp  = 

AK( 

CON 

cal 

IF( 

CM 
CAL 

r.AL 
on 

on 

STF 
STF 
CON 

On 

on 

STF 
STF 

con 
FnR 

RFT 

END 


(I  ,  J) 

TTNUF 
NDIM. 
?1  1  = 
=  5*1  + 

21  J  = 
J  +  l 
J)  =  AK 
NDIM. 
1?  1  = 
=  4*1  + 

22  J  = 
J  +  l 
J)  =  AK 
TINUF 
L  INV 
KFR.F 
L  MUL 
L  MAT 
L  MUl 
50  1  = 
50  J  = 
MOD<  I 
mod  (I 
TINUF 
60  1  = 
6C  J  = 
MnO(  I 
MOD(  J 
TINUF 
MAT  ( 
URN 


=TKT(I,J)+RB(K,I  )*AKT(K, J) 

EC. 6)    GO    TO    25 

1,5 

1 

1ST, 35 

<  JP) 

EQ.5 )    GO    Tn    25 

1,5 

1 

1ST, 35 

(JP) 

ERT     (NHIM, 1# 00-12. TKT,TKTINV,KER, NO  I M  ) 

0.21    WRITF (6, 1000) 

T(AKT,T«TINV,TTT,24,NDIM,N0IM) 

T{ AKT,AKTTt?4,N0IM) 

T(TTT, AKTT, ST FR B , 2 4, ND I M , 24 ) 

l»24 

1,24 

,  J)=0.CD0 

fJ)=STIFFF( I,J)-STFRB(I ,J) 

It  24 
It  24 

, JI=(STFMOD( I, J)+STFMOD( J, II )*0.5D0 
,  I)  =  STFMOO( I,  J) 

/,23H  MATRIX  K22  IS  SINGULAR) 


SUBROUTINE  M ATT ( A, AT t N 1 , N2  ) 
MATT  IS  AN  ONLINE  MATRIX  TRANSPOSE 
IMPLICIT  REAL*8( A-H.O-Z) 
DIMENSION  MM  »N2I  ,AT(N2,N1) 
on  ICO  1=1, N2 
on  100  J=1,N1 
100  AT( I , J)=A( J, I) 
RFTURN 
END 


ROUTINE, 


2 

10 


11 
12 
13 


14 


15 
20 


SUBROUTINE  INVERT(N,EP,A,X,KER ,NACT ) 
IMPLICIT  REAL*8( A-HtO-Z  ) 
DIMENSION  MNACT,NACT) , X(NACT,NACT ) 
DO  1  I=1,N 
00  1  J=1,N 
X( I , J)=0.000 
00  2  K=1,N 
X(K,K  )  =  1.0D0 
DO  34  L=l ,N 
KP  =  0 
Z=0.000 
DO  12  K=L,N 

IF(Z-DABS(A(K,L>  )) 11,12,12 
Z=DABS( A(K,L  )) 
K  P  =  K 

CONTINUE 

IF(L-KP)13,2C,20 
DO  14  J=L,N 
Z=A(L, J) 

A(L, J)=A(KD,J) 
A(KP, J)=Z 
DO  15  J=1,N 
Z=X(L, J) 
X(L, J)=X(KP, J) 
X(KP, J)=Z 
IF(DABS(A(L,L) )-FP ) 50, 50,30 


75 


*0 
31 


3? 

33 

35 

36 
34 
AO 


A2 
A3 


50 
75 


IF(L-N)31,3A,3A 

LP1=L+1 

DO  36  K=LP1,N 

IF( A(K,L)  )32,36,32 

RATIO=A(K,L) /A(L,L) 

00  33  J=LP1,N 

A(K,J)=A(K,J)-RATIO*A(L,J) 

00  35  J  =  1,N 

X(K,J)=X(K,J)-RATIO*X(L,J) 

CONTINUF 

CONTINUF 

DO  A3  1  =  1, N 

II=N+1-T 

00  A^  J=1,N 
S=0.000 

IF( II-N)A1,A3,A3 

1  IP  1  =  1  I  +  1 

00  A2  K=IIP1,N 

S=S+A(  II  ,K)*X(K 

X(  II , J)  =  ( X(  I  I,  J 

KER  =  1 

GO  TO  75 

KEP  =  2 

CONTINUE 

RETURN 

FNO 


tJ) 

I  -SI /AC  I  I  ,1  I) 


SUBROUTINE  MULT ( A, 8, P , N , M, L) 
MULT  IS  AN  ONLINE  MATRIX  MULTIPLICATION 

IMPLICIT  REAL*8( A-H,0-Z) 

DIMENSION  MN,M)  ,B(M,L)  ,R(N,L) 

00  20  1=1, N 

DO  20  J=1,L 

R( I , J)=0.ODO 

00  10  K=1,M 

R( I, J)=R(  I  ,J 
10  CONTINUF 
20  CONTINUF 

RETURN 

END 


ROUTINE 


)+A(I,K)*B(K,J ) 


1 

1000 


2 

2000 


3 

3000 


AOOO 
A100 

5 
5000 

6 
6000 


SUBROUTINE  FLAGfNO,*,*) 

IMPLICIT  REAL*8  (A-H,0-Z) 

RFAL*^  Kl  K?  K  3  I 

C0MM0N/B1/E,H1,H2,P0IS,RA01,RAD?,L,PHI1,PHI2,ALPHA, 
4THETA,KASF 

GO  TO  (1,2, 3, A, 5, 6), NO 

WRITE  (6,1000) 

FORMAT  (/,5X,'DATA  INCONSISTENT,  ARC  LENGTH  LESS  THAN1 
«, •  DIFFFPENCE  IN  RAOI  I' ) 

RETURN  1 

WRITE  (6,2000) 

FORMAT  (/,5X,'0ATA  INCONSISTENT,  COSPHI  IS  CALCULATED* 
#,'  TO  BE  GREATER  THAN  ONF') 

RFTURN  1 

WRITE  (6,3000) 

FORMAT  (/,5X, 'ELEMENT  IS  SPECIFIED  A  CYLINDER') 

RETURN  2 

WRITE  (6, AOOO) 

WIDTH=1.0D0/THETA 

WRITE  (6,A10C)  L, WIDTH 

FORMAT  (/,5X, 'ELEMENT  IS  SPECIFIED  A  FLAT  PLATE') 

FORMAT  (/,8H  LENGTH=,D15.8,7H  W I DTH= ,01 5 .8 ) 

RETURN  2 

WRITF  (6,5000) 

FORMAT  (/,5X,' ELEMENT  IS  SPECIFIED  A  TRUNCATED  CONE') 

RETURN  2 

WRITE  (6,6000) 

FORMAT  (/,'  SIN  PHI1  OR  SIN  PHI  2  IS  EQUAL  TO  ZERO,', 
*'  RIGID  BODY  MATRIX  WILL  HAVF  SINGULAR  POINTS, ',/, 
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#•  ELEMENT  STIFFNESS  MATRIX  OOFS  NOT 

#»BOOY  MOTION,  RELOCATE  NODE  POINTS  • 
RETURN  2 
FND 


CONTAIN 


RIGID  • 


SUBROUTINE 

STRFSS  FGPMS 
IMPLICIT  R 
REAL* 8  Kl, 
DIMENSION 
#P(6,24) fCO 
DATA  CORNE 

INITIALIZE  Wl 
DO  50  I  =  1  t 
DO  50  J=l, 

50  Wl  ( I , J)=0. 

CALCULATE  STP 
DO  100  1=1 
X=CCRNFR(I 
Y=CORNFR< I 
CALL  PST(P 
II=6*(  1-1 ) 
DO  100  J=l 
JJ=II*J 
DO  100  K  =  l 
100  W1(JJ,K)=P 
CALL  "ULT( 

RESET  Wl  MATR 

DO  200  1=1 

DO  200  J=l 

200  Wl( I ,J)=C, 

CALCULATE  STR 
DH  300  1=1 
X=CCPNEP( I 
CALL  CEIK 
11=6* ( 1-1) 
DO  3C0  J=l 
JJ=II+J 
DO  300  K=l 
KK=I I+K 
300  W1(JJ,KK)= 
CALL  KULTI 
RETURN 
END 


S 

EL 
FA 
K2 
SS 
RN 
R/ 
M 
24 
24 
OD 
M 
,4 
,1 
,2 
,X 

,6 


TRESS(TINV,SS, SN) 

EMFNT  STRESS  AND  STRAIN  RELATIONSHIPS  . 

L*B  (A-H,0-Z) 

K  3  I 

(2^,24  >,SN(  24,  24),T  INV(  24,2  4),W1(  24,24)  , 
ER(4,2) ,DFL(6,6) 

O.ODO, 2*1 .ODD, 3*0.000, 2* 1« ODO/ 
ATRIX 


0 

N  RELATIONSHIP,  SN  IS  STRAIN  MATRIX 


) 
) 
,Y) 


) 

INV,SN,24,24, 24) 
H  ZERO 


t24 

(J,K 

Wl  ,T 

IX  T( 

t24 

»?4 

one 

ESS  RELATIONSHIP,  SS  IS  STRESS  MATRIX 

»4 

,1) 

DEL,X  ) 

,6 

,6 

DEL( J,K) 

MltSNtSSf 24,24, 241 


SUBROUTINE  MERGER  (TOTK) 
C   MERGE  ARRANGES  INDIVIDUAL  ELEMENT  STIFFNESS  MATRICES  ON 
C   ON  THE  UPPER  HALFRAND  OF  THE  TOTAL  STIFFNESS  MATRIX  AND 
C   PLACES  APPROPRIATE  BOUNDARY  VALUES  IN  THE  LHAD  VECTOR. 

IMPLICIT  REAL**  (A-H,0-Z> 

DIMENSION  AK(24,24) ,LM(4) 

C0MMGN/MBN/BL0AD(2  94) ,NBC( 50, 7 ) ,NC0N<36. 51 

C0MMON/PAPM/NEL,NELT,NDPT,NPBC,NBAND,NF0 

DIMENSION  TOTKI 2^4,54) 

DO  4C0  N=1,NELT 

NTK=(NCCN(N,5)-1  )*4+l 

CALL  R0DISK(NTK,AK,576) 

LM(1 )=NCDN(N,1 )*6-6 

LM(2)=KC0N(N,2)*6-6 

LM( 3)=NC0N(N,3)*6-6 

LM(4)=NC0N(N,4)*6-6 

nn  aoo  1=1,4 

DO  ^00  J=l,4 

nn  400  k=i,6 

DO  400  LL=l,6 
II  =  LM(  I)+K 

JJ  =  LM( J)+LL-LM(  I  >-K  +  l 
IF(JJ)  40C, 400,300 
300  CONTINUE 

KK=6*I-6+K 
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LL 
AK  IS 

TCTK 

TO 


400 


500 


60C 
700 
800 


CO 
DO 
IP 
DO 
IA 
IF 
DO 
TO 
TO 
3L 
II 
LL 

no 

LL 
II 

IF 

TO 
CO 
CO 
CO 
RE 
EN 


L  =  6 
EL 
IS 
TK( 

NTI 

fiO 
T=( 

70 
0=1 
(NB 

50 
TK( 
TK( 
OAC 
=  0 
U  =  N 

60 
L  =  L 
=  IA 
<  II 
TK( 
NTI 
NTI 
NTI 
TUR 
D 


*J-6+LL 

EMENT  STIFFNESS  MATRIX 

THF  TOTAL  STIFFNES  MATRIX 

II,JJ)=TOTK< I  I, JJ)  +  AK(KK,LLL) 

NUF 

0  I=1,NPBC 

NBC( I  ,1  )-l  >*6 

0  J=l,6 

PT  +  J 

C( I , J+l  ).E0.1 )  GO  TO  700 

0  K=1,NBAND 

I AC,K)=0.0D0 

IAD,  l)  =  l.ODO 

(  IAD)=0.0D0 


BAN 
0  I 
LM 
D-L 
.Lp 
lit 
NUE 
NU^ 
NUF 
N 


D-! 
L=1,LLU 


.0)  GO  TO  600 
LLL)=0.0D0 


SUBROUTINF  8ANS0L( A,B,NN.MM) 
C   BANSOL  SOLVES  THE  TOTAL  FORCE-DISPLACEMENT 

IMPLICIT  RFAL*R  (A-H,0-Z) 
C         NN  IS  THE  NUMBER  OF  EQUATIONS 
C  mm  IS  THE  HALF  BAND  WIDTH 

DIMENSION  A(2<?4,54),B(294) ,C(54) 

N  =  0 
100  N=N+1 

B(N)=B(N)/A(N,1) 

IF(N-NN)  150,300,150 
150  DO  200  K=2,MM 

C(K)=A(N,K) 
200  A(N,K)=A(N,K)/A(N,1) 

DO  260  L=2,MM 

I=N+L-1 

IF(NN-I)  260,240,240 
240  J=0 

DO  250  K=L,MM 

J  =  J-M 
250  A< I,J)=A( I,J)-C(L)*A(N,K) 

B(I)=B(I )-C(L)*B(N> 
260  CONTINUF 

GO  TO  100 
?00  N=N-1 

IF(N)  350,500,350 
350  DO  400  K=2,MM 

L=N+K-1 

IF(NN-L)  400,370,370 
370  B(N)=B(N)-A<N,K)*B(L ) 
400  CONTINUE 

GO  TO  3C0 
500  RETURN 

END 


FQUATIONS. 


SUBROUTINE  PLOAD( P , PVEC ) 
IMPLICIT  REAL*8  <A-H,0-Z) 
REAL**  L 

DIMENSION  PVEC (24) ,QVEC( 24),XG( 8  ) ,  WT  (  8  ) ,C0EF(4) ,F(4) 
#TINV(24,24) ,TINVT( 24,24) 

C0MM0N/B1/E,H1,H2,P0IS,RAD1,RAD2,L,PHI1 , PHI  2, ALPHA, 
#THETA,KASE 
INTEGRATF  TERMS  OF  THE  PLOAC  VECTOR, 
ZEROS  OF  GAUSSIAN  8  POINT  FORMULAS  ARE 
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XG 


XX1  =  . 

XX2  =  . 

XX3=. 

XX4  =  . 
COEFFICI 

AA1=. 

AA2=. 

AA?  =  . 

AA4=. 

CONTA 
COFFFICI 

WT(1  ) 

WT(2) 

WT(3) 

WT(4) 

WTC5) 

WT<6> 
WX(7I 

A=0.5 
XG<1> 
XG(2) 
XG(3) 
XG<4) 
XGCS) 
XG(6) 
XG(7) 
XG  (  8 ) 
DG  10 
F(  I  )  = 
DO  20 
QVFC( 
X  =  XG( 
WTW  =  W 
F  (  1  )  = 
F(  2)  = 
T(3)  = 
F<4)  = 
CONTI 
CALCULAT 

on  30 

300  COFF( 
LOAD  T[;p 
DO  40 
DO  40 
K=<  J  + 
QVFC( 
CONTI 
TRANSFOR 
CALL 
CALL 
RETUP 
END 


100 


200 


400 


9602S 
79666 
52553 
18343 
FNTS 
10122 
22238 
3137C 
?6?t£ 
IMS  A 
FNTS 
=  AA1* 
=  AA2* 
=  AA3* 
=  AA4* 
=  WT(4 
=  WT(  3 
=  WT(2 
=  H(1 
DO 

=  A  +  XX 
=  A  +  XX 
=  A  +  XX 
=  A  +  XX 
=  A-XX 
=  A-XX 
=  A-XX 
=  A-XX 
0  1  =  1 
O.ODO 

0  1  =  1 

1  NO. 
I  ) 
T(  I  ) 

(1)  + 
<?)  + 
C)  + 
(4)  + 

NU^ 

E  COE 

0  1=1 

1  )=P* 
MS  OF 
0  1  =  1 

0  J=l 

1  )*4  + 
K)=CO 
NUE 

M  OLC 
TFORM 
MULT( 
N 


■9753623D0 
1362674D0 
163289900 

■  956493000 
POS  APE 
90376?6D0 

,S337447no 

;7788729D0 
7836193D0 

i«=  P  ZEROS 


564« 

774' 

C90 

424« 
7FI 

362 

844' 

458 

633 
TH<=  8  ZEROS  AND 
ALL  ARE  SHIFTED 

500 

5D0 

5D0 

^DO 


WT  CONTAINS  ALL  THE 
TO  THE  RANGE  0  TO  t 


93 
64 
->4 
46 
OF 
85 
10 
66 
37 
LL 

0. 
C. 
C. 
0. 

) 
) 
) 
) 

1*0. 5D0 
2*C.5D0 
3+0. ^DO 
4*0.500 
4*C.?00 
3*0. ^DO 
2*0, 5D0 
i*C.5D0 
t4 

t8 

ODO 


WTW*R(X)*X**3 
WTW*<? ( X)*X**2 
WTW*R ( X)*X 
WTW*R (X) 

FFICIENTS    FOR     INTERGRALS 

,4 

L*ALPHA/(5.0D0-I ) 

OLOAD  VECTHR 
,4 
t4 
I 
EF (I >*F< J) 

AD    VECTOR    TO    GET    PLOAD    VECTOR 
(TINV,TINVT) 
TINVT,QVFC,PVFCf 24,24,1 ) 


SUBROUTINE 
C  WRDISK/PDriSK 
C       DISK. 

INTEGER    LAST/0/ 

REAL*8    NAME<2)/,WRDISK« 
C       REMOVE    THE    FOLLCWING    CAPO 

RFAL*8     A 

DIMENSION    A(l) 

DEFINE  FILE  7( 
C   REPLACE  368  'WITH 

IFfNTRACK  .GT. 

IF(NTPACK  .LT. 

N  =  NTPACK 
C   THF  MAXIMUM  NUMBER 

NRL=13 

N=N*NRL+1 

IF  (NCT.GT.  46)G0  TO  50 


WRDISMNTPACK,  A,NCT) 
READS  AND  WRITES  ELEMENT 


INFORMATION  ON/OFF 


WHEN  USING  SINGLE  PRECISION 


1500, 368, E, I ) 

184  FOR  SINGLE  PRECISION 

1999)  GO  TO  900 

0)     NTRACK    =<LAST+9)/10 

OF    WORDS     IN    A    OR    B    MUST     BE    .LF.     NRL*46 
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WRITF  (7«N,1000)  (  A(  J )  , J  =  l ,NCT ) 

IF  (LAST  .LT.  I )  LAST  =  I 

IFUAST.GT.  19999)LAST=0 

RFTURN 
50  JI  =  47 

WRITF  C«N,1000)  ( A( J)  ,J  =  1  ,46) 
75  JF  =  J  I  ♦  45 

IF  (  JF  .GF.  NCT)  GO  TO  100 

WRITE  (7M,1000)  ( A( J)  ,  J  =  JI , JE) 

JI  =  JI+46 

GO  TO  75 
100  WRITE  (7M.1000)  (  A  (  J  >  ,  J  =  J  I  ,  NOT  ) 

IF(  I.GT.LAST)  LAST=I 

IF(LAST.GT,     19Q99)LAST=0 

RCTUPN 

ENTRY  ROOISK(NTRACK,B,NCT) 
C   REMOVE  THE  FOLLCWING  CARC  WHEN  USING  SINGLE  PRECISION 

REAL*S  B 

DIMENSION'  B(l) 

NRL  =  13 

ITINTRACK  .GT.  1999)  GO  TO  910 

N=NTRACK 

N=N*NPL+1 

IF  (NCT  .GT.46)  GO  TO  150 

READ  <7»N,1000)  (B(J),  J=1,NCT) 

RETURN 
150  READ  (7'N,1000)  ( B ( J ) , J=l , 46 ) 

JI  =  47 
175  JF  =  JI  +  45 

IF  (JE  .GE.  NCT)  GO  TO  200 

READ   (7*1,1000)  (B( J),J  =  J  I,JE  ) 

JI  =  JI  ♦  46 

GO  TO  175 
?00  READ   (7M,1000)  ( B( J ) , J  =  J  I , NCT  ) 

RETURN 
1000  FORMAT  (46A8) 
C   PEPLACF  FORMAT  WITH  (46A4)  FOR  SINGLE  PRECISION 
900  K=l 

GO  TO  <5?C 
^10  K=2 

920  WRITE(6, 1920)NAME( K) ,NTRACK 
1920  FORMATCO  ERROP  IN  CALL  OF  SA6,',  NTRACK  TOO  LARGE, (M 
1THAN  200C)'  /  •   NTRACK  =»,U0) 

STOP 

END 


100 
200 


SUBROUTINE  FCRML 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  K1,K2,K3,L 

DIMENSION  RVEC(24) ,LM(4) 

C0MM0N/MBN/BL0AD(2  94) ,NBC( 50 , 7 ) , NCON ( 3 6 , 5 ) 

C0MM0N/PARM/NEL,NFLT,N0PT,NP3C, NBAND,NEQ 

00  200  I=1,NELT 

NTK=NCON( I,5)*4 

CALL  RDPISK(NTK,PVEC,24) 

LM( 1)=6*NC0N(I ,l)-6 

LM(2)=6*NC0N( I,2)-6 

LM(3)=6*NC0N(I ,3  )-6 

LM(4)=6*NC0N(I,4)-6 

DO  100  K=l,4 

II=LM(K) 

IP=6*(K-1) 

DO  100  J=l,6 

IB*II+J 

IIP=IP+J 

BLOAD(  IB)  =  BLCAD(  IB  >+PVEC( IIP) 

CONTINUE 

RETURN 

END 
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SUBROUTINF  PSSN 
C   THIS  SUBROUTINE  RECOVERS  ELEMENT  STRESS  kfifi    STRAIN  VALUES 
C   FROM  DISK  AND  ARRANGES  THEM  IN  RESULTANT   STRESS  AND 
C   STRAIN  MATRICES  FDR  THE  NCOF  POINTS. 

IMPLICIT  REAL*8  (A-H,0-Z) 

CCMM0N/MRN/8L0AD(294)  ,NBC<^0,7)  ,NCnNH6,5) 

COMMON /PARM/NEL,NELT, NO PT,NP8C , NBA  NO, NFQ 

COMMON/ SS  SM/SS ( 24,24) ,SN( 24,24) 

DIMENSION  UEL(24)  , SNL ( 24 ) , SSL ( 24 ) , R ESN ( 49,  "M  , 
#RCSS(49,6)  ,LM(4) 
C   INITIALIZE  RESULTANT  STRAIN  AND  RESULTANT  STRESS  MATRICES 

DO  100  1=1,49 

DO  100  J=l,6 

RESN( I  ,  J)=0.000 
100  «ESS(  I  ,J)=0. COO 

DO  110  1=1, 49 
110  RESN(  I  ,7)=0.000 

00  400  I=l,NFLT 

NTK=4*NC0N( I ,5)-2 

CALL  ROOISK  (NTK,SN,576) 

NTK=NTK+1 

CALL  RDOISK  (NTK,SS,576) 

LM(  1)=6*NC0N(I  ,  1  )-6 

LM( 2)=6*NC0N( I , 2)-6 

LM(3)=f*NC0NUf3)-6 

LM(4)=6*NCCN(I  ,4)-6 

DO  200  J=l,4 

JAD=LM(J) 

IUEL=( J-l) *6 

DO  200  K=l,6 

KJA=JAD+K 

IU=IUEL+K 
200  UEL( IU)=BLOAD(KJA) 

00  210  111=1,24 

SNLUII  )=C.0D0 

00  210  JJJ=1,24 
210  SNLIII  I  )  =  SNL(  II  I  )  +  SN( I  II,JJJ)*UEL< JJJ) 

00  220  111=1,24 

SSL (III  )  =  C.0D0 

DO  220  JJJ=1,24 
220  SSL( III )=SSL(III )  +  SS(III , J  J  J > *UEL( J J J) 

DO  ^00  11=1,4 

KJA=NCCM  I,  II) 

IUEL  =  (  II  — 1 )*6 

00  250  JJ=1,6 

IU=IUEL+JJ 

RPSNIKJA,  JJ)  =RESN(KJA,  JJX-SNLI  IU) 
250  RESS(KJA,JJ)=RESS(KJA,JJ)*SSL(IU) 
300  RFSN(KJA,7)=RESN(KJA,7)+1.000 
400  CONTINUE 

DO  500  I=1,NDPT 

00  500  J=l,6 

RESN( I , J)=RESN< I ,J)/RESN(I ,7) 
500  RESS( I ,J )=RESS(  I, J  )/PESN(I ,7) 

WRITE  (6,1000) 
1000  FORMAT ( 1H1 ,/// , 27X, • STRAINS' ) 

WRITF  (6,1500) 
1500  FORM AT <//,?7x,» JOINT* , 5X , • SNX' , 9X, • SNY • , ^X , • SNXY • f SX, 
#«SMX« ,9X,  'SMY' ,9X,  'SMXY*  ) 

DO  600  I=1,NCPT 
60C  WRITF  (6,2C00)  I , ( RE SN ( I  , J ) , J  =  l  ,6 ) 
2000  F0RMAT(/,19X,I 12,1P6012.4) 

WRITE  (6,2500) 
2500  FORMAT! 1H1 f /// f 2 7X,« STRESSES1 ) 

WRITE     (6,3000) 
3000    FORMAT (//,?7X,«  JOINT'  ,6X, 'NX' , 10X,  «NY»  , 
IH9Xt  *HX*t  lOXt  •MY*  »10X,  *MXY*  ) 

00    700     I=1,NDPT 


700  WRITF  (6,2000) 
RETURN 
END 


I,(RESS(I,J), J=l ,6) 


1000 

I5cc 


100 
2000 


SUBROUTINE  OUTPUT 
OUTPUT  WRITES  NODAL  DISPLACEMENTS. 
IMPLICIT  REAL**  (A-H,0-Z) 
RFAI  *fl  Kl  K?  K "^  I 

COMMON /MBN/BLOAD( 2 9M  ,NBC( 50 , 7 ) , NC ON ( 36 , 5) 
COMMON /PARM/NEL,NELT, NO PT,NPBC,  NBAND,NEQ 
WRITE  (6,1000) 

F0RMAT(1H1,///,27X, •DISPLACEMENTS* I 
WRITF  (6,1500) 

FORMAT ( //,?7X,»J0INT«6X,«U»,11X 
#9X,«W,Y«  ,QX,  'W^Y'  ) 
no  100  1=1 ,NCPT 
IL=6*(  1-1)4-1 
IU=IL45 

WRITE  (6,2000)  I , ( BLOAD( J )  , J= IL ,  IU ) 
FORMAT (/,19X, I 12,1P6D12.A) 
RETURN 
END 


,'V«  ,11X 


•  ,10X,  «w 
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